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ABSTRACT 

An improved mathematical model using the Vrentas-Duda theory 
for diffusion coefficients is developed for solid state 
polymerization (SSP) . This model is applied to nylon 6, and 
best-fit values of the parameters are obtained for this system 
using the Box-complex analysis on available experimental data. The 
sensitivity of results to variations of these parameters is 
studied. It is found that the ring opening reaction can easily be 
omitted, and the number of parameters thereby decreased. Detailed 
quantitative results are obtained to study the effects of changing 
the important operating conditions on SSP, e.g., intermediate 
remelting of nylon 6 powder, value of the water concentration (or 
level of vacuum) in the vapor phase, leaching of monomer and water 
before SSP, size and degree of crystallinity of polymer particles, 
etc. The diffusional effects on the rate constants are found to 
depend upon a complex interplay of several important factors like 
the generation of monomer and water by reaction, removal of water 
by macro-level diffusion, changes in the intrinsic rate constants 
due to decrease in the concentration of polymer molecules, etc. 



CHAPTER 1 

INTRODUCTION 


1 


It is often difficult to obtain extremely high molecular 
weight polymers by conventional melt polycondensations because of 
high melt viscosities and limitations posed by equilibrium. 
However, it is possible to raise the degree of polymerization, DP, 
of several step growth polymers by heating them in vacuum or in an 
inert gas at temperatures below their melting points but above 
their glass transition temperatures. This method, known as solid 
state polymerization (SSP) is quite commonly used in industry. 
Applications involving stringent operating conditions (e.g., tire 
cords) and certain processing techniques like blow molding and 
extrusion require such polymers. Both batch as well as continuous 
operations are being practiced in industry for this purpose, 
involving the use of fixed bed, rotary drum and fluidized bed 
reactors. 

Studies on SSP go back to the early forties. Most of the 
early information is in the form of patents. Work on SSP started 

getting reported in the open literature about two decades ago. 1-4 

5 

Gaymans et al. carried out a systematic experimental study on the 
SSP of nylon 6 , while Griskey and Lee 4 and Chen et al . 1 provided 
some experimental data on nylon 66 and on nylon 6-10 and 
polyethylene terephthalate respectively, under a variety of 
conditions. Pilati and Fakirov have provided excellent reviews 
of the work in this area. 

Most of the early studies provided semi-empirical models to 
explain experimental data. Kumar et al. proposed a 
phenomenological model having a molecular basis to explain some 
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important trends observed in the SSP of condensation polymers. 

9 , 

They adapted the model of Chiu et al. for the gel effect in free 

radical polymerizations, to SSP. More recently, Kaushik and 

Gupta 10 extended and improved this model and applied it to the SSP 

of nylon 6. They obtained curve-fitted values of the parameters in 

. 5 

their model, using the experimental data of Gaymans et al. In the 
present study, we present a new model, which uses the free volume 
theory of Vrentas and Duda 11 to account for the changes in the 
diffusion coefficients more accurately, using parameters which can 
be obtained directly from independent experiments which do not 
involve reactions. This model, applicable to any reversible step 
growth polymerization in the presence of diffusional limitations, 
is applied to the specific case of the solid state polymerization 
of nylon 6 in this work. An optimal parameter estimation package 

is again used to obtain the best-fit values of the parameters used 

• • • 5 

m this model, using the experimental data of Gaymans et al. 
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CHAPTER 2 

FORMULATION 

There are three major reversible reactions in the kinetic 
scheme of nylon 6 polymerization : ring opening, polycondensation 
and polyaddition. These are shown in Table I. The reaction most 
susceptible to diffusional limitations is the forward step of the 
polycondensation reaction, since it involves the diffusion of two 
large molecules towards each other. In SSP however, it is expected 
that the diffusion of smaller molecules like water and 
e-caprolactam can also be affected. The present model takes this 
also into consideration. Thus, the reactions which are considered 
to be diffusion controlled in the present study are the forward 
and reverse steps of polycondensation, as well as the forward step 
of the polyaddition reaction, since at least one large molecule is 
involved in these three reactions. The backward steps of ring 
opening and polyaddition reactions are not diffusion controlled, 
since they involve only a single species. Similarly, the forward 
step of the ring opening reaction is assumed not to be diffusion 
controlled since it involves the diffusion of two small molecules. 

We now model the variation of the forward rate constant of 
the polycondensation reaction associated with diffusional 
limitations. The approach used is similar to that suggested by 
Achillas and Kiparissides. 13 ' 14 Fig. 1 shows a polymer chain, P^, 
with its functional A group at the center of a hypothetical sphere 
of radius r^ Within this sphere, the effective concentration 

of polymer molecules is [P] m . At a large distance, r Q , from the 
central functional group, the concentration of polymer molecules 
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TABLE I 

KINETIC SCHEME FOR POLYMERIZATION OF NYLON 6 


Ring Opening 


c ± + w 



Polycondensation 


P n + 


m 



P 


1 


n+m 


+ w 


Polyaddition 


P n + 



P 


n+1 






6 


approaches the (local) bulk value, [ P] . The reaction between the 
central functional A group and a B group of another polymer 
molecule, P , far away, occurs only after the latter comes in the 
vicinity of A by diffusion. The diffusion rate can be easily 

Q 

shown to be independent of the radial location, and is given by 


4 


2 


it r 


d[P] 
d r 


CD 


where is the micro-level diffusivity of the polymer molecule in 
the medium. The following boundary conditions apply 

at r=r m,l ‘ t p ) m < a > 

at r = r D — > » [P] = [P] fa (b) (2) 

Eqns. (1) and (2) give 


4 71 D p r m,l < ^b-tSV 


= C. 


(3) 


The rate of consumption of functional group B within the sphere of 
radius r . because of reaction between two large molecules 

Id, 1 

(forward reaction of polycondensation step) is 

4- n 4,1 < k 2,o [ p ]» [ p ) b > - C 1 < 4 > 

r 

where k 2 Q is the intrinsic rate constant for this reaction. Eqns. 
(3) and (4) give the expression for the effective polymer 
concentration at r = r . as 

HI, 1 
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. 2 

The overall rate of this reaction can also be written as ^[PJ^, 
where k 2 is the apparent rate constant for this reaction, 
accounting for diffusional limitations as well. The two 
expressions for the rate give 


k 2 C p lb “ k 2,o f p lm t p lb <«> 

On substituting the expression for [P] from Eqn. (5) in Eqn. (6) , 
we obtain 




( r m,l 
3 °P 


) [ p 3 b 


(7) 


where [P]^, the local bulk polymer concentration, can be rewritten 
in terms of the molecular weight distribution and its moments 
(see nomenclature) as 


C p 3 b 


£ 

n=l 


tVb - 


o,b 


( 8 ) 


An analogous mathematical treatment for the backward step of 
the polycondensation reaction and the forward step of the 
polyaddition reaction leads to 


*2 


TT 7- + < ) C p 3 b 


2 , o 


(a) 


w 


TET- + < TTT > tP3 b 

3,o m 


(b) 


(9) 


where D and D are the effective micro-level dif fusivities of 
w in 

water and monomer in the medium. 

The free volume theory of Vrentas and Duda 11 can now be used 
for the dif fusivities. The effective diffusivity, D^, of polymer 
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molecules can be written as 


D po “n 2 exp(_ R^T > exp < - * ) 


( 10 ) 


where 


A * A £ 

p <p V p <p V 

K nr m m rc? 


Z = 


•13 


s - -g. .. + p <p V* 

? 23 p p ”p p 


} 


A* 

m r m V m * fm 


P* V_ V 


+ p <p V* V„ + p <p V* V, 
K s r s s fs K p r p p fp 


( 11 ) 


The same notation has been used here as in Ref. 13 (also see 
nomenclature) . A semiempirical procedure is now used to simplify 
Eqn. (10) further. A reference state is defined as <p > 0, <p » 

Xu S 

0, and > 1. Eqn. (10) can be written in terms of the reference 

state (ref) as 


1 

k. 


"2,o 


'a, 1 


3 D. 


Pr ref J 


(D P 7 D P,ref> 


" 2,0 


~m, 1 


3D p, ref 4 n,ref' 


M n * 0 ex P ( X 


~ * re f ) 


where 


'ref 


"2,o 

J_ 

V.* 

fp 


+ ©t (T) M n X Q exp ] 


ref 


( 12 ) 


(13) 


13 

The parameter, (T) , is a function of temperature alone for a 
given sample (it depends on the degree of crystallinity through 
the parameter ref ) • In a similar manner, Eq. (9) can be 


written as 
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+ 0 2 (T) exp [€ 23 (Z “ * ref ) ] ( a ) 

+ ® 3 (T) X q exp t€ 13 (X " X ref )l < b > (14) 

Table II summarizes the final equations for the apparent rate 
constants. Numerical values of and correlations for the various 
parameters required have been compiled from various sources and 
are given in Table III along with references. ' These 

parameters have been determined experimentally using independent 
measurements on nonreacting nylon 6 systems. The glass 

transition temperature of nylon 6 has been assumed to be a 
function of its spatially averaged value of the number average 
chain length, ji^, as assumed earlier 8 ' 10 . 

The micro-level diffusivity of water leads to a difference in 

the water concentrations at r = r , (where the concentration is 

m, 2 

[W] m ) and r = r D » » (where the concentration is [W] b ) . This 

gradient is accounted for in the expression for the apparent rate 

constant, k', in Eqn. (14a). It may be pointed out that the bulk 

concentration of water, [W] b , itself is not constant, but varies 
with location due to the macro-level diffusion of water through 
the polymer pellet to the flowing inert gas outside. The 
diffusional resistance of W in the gas phase is neglected in this 
formulation. However, it is expected that there is considerable 
resistance for the diffusion of water within the pellet. We 
consider the reacting pellets to be spherical with the 

condensation byproduct, W, diffusing radially outwards. The 

macro-level diffusivity, D w , 



of W through the 'solid' nylon 6 



TABLE II 

APPARENT REACTION RATE CONSTANTS 


1_ 

k. 


+ 0, (T) u: A 


2 , o 


n o 




exp - * + x 


ref 


v ' 
K 2 


w 

K 2 , o 


+ 0 2 (T) 


exp[? 23 ( - X + X ref )_ 


1_ 

k. 


±r* e 3 < T > 

3,0 


exp[e i3 ( - * + x ret )_ 


v • k ' 
K l,o' K 1 


= k' = k- /K 

1,0 1 , 0 ' 1 


V ' 
K 3 


kf ^ = k /K 0 . k' = k 0 /K n 
3,0 3,0 3 , 2,0 2,0 2 


'ref 


■13 


{ 


A* 

P**™ v. 


A * 

P e 


m T m 'm , r 's 7 's's , , t A * ] 

+ + Pr, | 


■13 


•23 


P P P 


P m <P m V* V- + p c <*> V* V- 
m m m fm 's s s fs 

r 


V- 

fp 

A * 

V (MW ) 
m v irr 

* * 

V M. 

P DP 


P_ V* V- 
P P P fP 


*23 


V* (MW g ) 
** 

V M . 

P DP 


m 


1 / v m < T > ; Pp - 


l/Vp“ ( T); p s 


1 / V S “(T) 


Contd. . . 
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(Table II contd...) 


A 

[C. ] (MW ) V° 
L 1 m m 




A 



A 




A 

[Cl] 

<“V 

V m 

m 

+ x 1 

(MW m ) 

v° 

p 

4- 

[W] 

(MWg) 






A 









X, (MW ) V 
1 ' m' 

O 

p 







A 



A 





t C l] 


v° 

m 

+ 

(MW m ) 

v° 

p 

+ 

[W] 

(MW s ) 

v° 

S 






A 








[W] 

(MW s ) 

V s 







A 



A 




A 

[C,] 

(MW m ) 

v° 

m 

+ A 1 

<"V 

V P 

+ 

[W] 

(MW s ) 



p : polymer 

s : solvent (water) 

m : monomer (e-caprolactam) 

local bulk values of concentrations and moments used for all 
cases (subscript b not indicated) 




TABLE III 

PARAMETERS USED TO CALCULATE THE APPARENT REACTION RATE CONSTANTS 


<“V 

= 

0.11316 kg/mol 


(MW s ) 

- 

0.018 kg/mol 


A * 

V 

= 

8.64X10 -4 m 3 /kg 

(Ref. 16) + 

m 




V m °(190°C) 

= 

11.325X10 -4 m 3 /kg 

(Ref. 16) + 

A * 

V 

D 

= 

8.123X10* 4 m 3 /kg 

(Ref. 17) + 

A Q 

V (T) 

P ' 

= 

V [1 + v f P <T)] 


* * 

V 

s 

= 

9.69xl0 -4 m 3 /kg (Water) 

(Ref. 15) + 

A o o 

V (190 C) 

= 

11 . 429xl0~ 4 m 3 /kg (Water) 

(Ref. 18) 

V fm ( 19 0°C) 

= 

0.31033 

(Ref. 16) + 

V fs (190°C) 

- 

0.1795 

(Ref. 15, 18 ) " 

V fp (T) 

SB 

0.025 + (a x - « g ) (T - T g ) 

(Ref. 19, 21) 

MW 


2 


m 

M. 



(Ref. 20 ) 

DP 


3 


r 

= 

1 

(Ref. 13, 14) 

1 


, K 


T 

= 

1 | / s v 1 

T K 2 } - 

(Ref. 8,19) 

g 


V u n 


V 

= 

350 K 

(Ref. 17) 



-3 -1 


a l 

= 

10 K 

(Ref. 17, 21) 



-4 -1 


a 

g 

== 

4.5 X 10 K 

(Ref. 17, 21) 

K S / T ga> 

= 

0.03 K -1 

(Ref. 10) 

K 1 

SB 

exp [ (-32.989 - 8024 . 9/T) /IR '] 

(Ref. 22) 

K 2 

= 

exp [ (3.9846 + 2 . 4 877xl0 4 /T) / IR ] 

(Ref. 22) 

K 3 

= 

exp [ (-29.06 + 1.6919xl0 4 /T)/IR ] 

(Ref. 22) 

OR 

= 

8.314 J/mol-K (Gas Constant) 



+ estimated values 



sphere is assumed to be constant. The mass balance equation for 
the spherical pellet with diffusion of W in the radial direction 
can be written as 


m 

a t 


^w r 2 a r y a 


m 


JR 


W 


( 15 ) 


where !R is the rate of 'generation' of water by chemical reaction 
w 

(ring opening and polycondensation reactions) and is written as 
R w = “ k i CC 1 ] tW] + k i [P l ] + k 2 X o “ k 2 (A 1 “ A o> [W] (16) 


The complete set of mass balance and moment equations are 
given in Table IV. The initial conditions can be taken as a 
uniform distribution of W in the pellet, i. e., [W] = [W] Q for 
all values of r at t = 0. The boundary conditions are easily 
written and are incorporated in Table IV. 

The equations in Table IV can be integrated along with the 
equations for the apparent rate constants (Table II) for a given 
set of parameter values, using the finite difference technique (11 
grid points) in the radial direction, and Gear's technique. The 
intrinsic rate constants, , are assumed to be of the following 
form 10 


I/O 


" <i 


i~l / 2 , 3 


( 17 ) 


where £ 2 and are the parameters (functions of temperature) 

which are curve fitted in this study. The dependence of k. on 

1 f o 

the concentration, X Q , of the acid end groups at any location is 

. . 12 

assumed for the SSP of nylon 6 because it is well known that 
these groups exert a catalytic effect in melt polymerization of 


TABLE IV 


FINAL EQUATIONS FOR SOLID STATE POLYMERIZATION OF NYLON 6 


1 . 


2 . 


k l C 1 W + k l P 1 ~ k -> A ~ + k ^ (A ~ “ Pl) 


1 o 


ffl „ 

at 

8P- 

at - k l C 1 B - k l P 1 - 2k 2 P 1 X o + 2k 2 W < X o - P l> 


k 3 P 1 C 1 + k 3 P 2 


4. 


° - k. C. W - k^ P. - k_ A 2 + k_ W( A. - A ) 
11 11 2 o 21 o 


= t C, W - k, P, + k, C, k - k_ (A - P n ) 
11 11 310 3 0 1 


3A 

at 

dx 

at 

^A_ t _ 1 J 

5. at = k i c i w “ k i p i + 2k 2 A 1 + I k 2 w (A i " A 3 ) 
+ k 3 C l( X 0 + 2A 1 ) + kj (A 0 - 2 X 3 + Pj) 

6 - “ ' k i c i w - k i p i + k 2 x o - k 2 W < X 1 ' X o> 


+ D [ 1“ + 2 ] 

w L Q 2 r 3r J 
or 

Initial and Boundary Conditions 

7. t = 0, values from Table V, independent of radial 


r = 0, 

r = R, 

Closure Conditions 


position 

aw 


ar = 0 

w * w 


8 

9. 


P 2 P 1 

A 3 = A 2 (2A 2 A o - A 2 )/ (A 1 A o ) 


* Brackets, [ ], for concentrations (local bulk values) omitted 
for brevity 


c -caprolactam, and a similar effect could be present in SSP as 
well. However, we have also carried out optimal parameter 
estimation studies assuming constant values of k. : 

± t o 

k. = ; i = 1/2,3 (18) 

1,0 1 

The equilibrium constants for the three reactions in Table I 

. . 12 

are assumed to be the same as those for melt polymerization. 
Thus, the intrinsic rate constants for the reverse reactions are 
given by 


V > 

_ k i/C 

; i - l, 2, 3 


(19) 

1 , O 

K i 


Since 

there are 

no diffusional limitations for 

the 

forward and 

reverse steps of 

the ring opening reaction, as 

well 

as for the 

backward step of 

the polyaddition reaction, we have 


k l 

= k l,o 

(a) 



k i 

- k i:o 

(b) 



k 3 

- k 3: 0 

(c) 


(20) 


The parameters © 1 , © 2 and 6 3 have been assumed as curve-fit 

parameters in this study (along with i=l,2,3). It may be added 

. . . . . 13 14 

that Achillas and Kipanssides ' have suggested the use of the 
theory of excess chain end mobility ' for estimating these 
parameters. However, they were forced to introduce a curve fit 
parameter, j , in this theory to obtain good agreement with 

C / O “ 

experimental data on the isothermal polymerization of methyl 
methacrylate (MMA) . Since a curve-fit parameter is necessarily to 
be used, we decided to use it for the entire group of variables 



clubbed together in the 9^ s. A similar procedure was found to be 

2 6 

successful in our recent work on the modeling of MMA 

polymerization in semi -batch reactors. 

The macro-level diffusivity, D^, of water through the reaction 

mass, is an important physical parameter in this work. Strictly 

speaking, tt> w would be expected to decrease with the degree of 

crystallinity, and so would be a function of the reaction time, 

since the polymer crystallinity increases as its molecular weight 

increases. Also, the nonuniform molecular weight profile within 

the polymerizing particle would imply that the crystallinity 

would be a function of the distance from the surface of the 

particle. Thus, the diffusivity is not only time-dependent, but 

also position-dependent. In this study, however, we have assumed 

D w to be constant at an average value. From the available data 27 

on the diffusivity of water through nylon 6 at a few temperatures, 

we estimated D w at the reaction temperature (190 C) to be 10 
2 

m /hr. Use of this value does not give a good match with 

experimental results. This value is much too high. On closer 

scrutiny, it was found that the sample (film) of nylon 6 used to 

generate the correlation for E> w had a crystallinity of around 
28 5 

30%. Gaymans et al. have not provided any information on the 

crystallinity of the samples they used for generating 

experimental results on the SSP of nylon 6. We guess that the 

crystallinity of their samples (molecular weights around 20,000 - 

30,000) could be as high as 75%. This would imply that D w would be 

•"8 2 

much lower than the value of 10 m /hr predicted using the 
correlation in Ref. 27. Because of this uncertainty, we decided to 
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treat D w as a parameter in this simulation, alongwith and 0^ 
(i=l,2,3) . 

In order to obtain the initial conditions to simulate the SSP 

5 . 

at 463 K for which some experimental data are available m the 

open literature, we first simulated the previous stage, namely, 

. . 12 

melt polymerization of e-caprolactam using [W] = 0.16 mol/kg, 

[C 1 ] 0 = 8.8 mol/kg, T=513 K, until a desired value of the number 

average chain length, h n , was obtained. This provided values for 

all the concentrations and moments at the end of the liquid phase 

polymerization. These values were taken to be the starting values 

(at t=0) for the SSP, except for the concentrations of 

e -caprolactam and water. [C^] and [W] at t=0 for the SSP, were 

taken as zero to duplicate their extraction out of the pellets 
. 5 

produced experimentally by Gaymans et al. from melt polymerized 
nylon 6. The starting concentrations used for simulation of SSP of 

nylon 6 are listed in Table V for their different experimental 

5 , ... 

runs (corresponding to differnt values of the initial number 

average chain length, Q ) • It must be emphasized that the 

, 5 

complete details on the conditions used by Gaymans et al. for 

producing the initial samples for SSP are not available in their 

paper, and this is why we have tried to generate initial values 

using the above technique, trying to match the only variable 

(h n Q ) whose value is mentioned. 

An optimal parameter estimation computer program has been 

used to obtain the best-fit values of the seven parameters (£^, 

0^; i=l,2,3; E> w ) using the isothermal experimental data on nylon 6 

29 30 

SSP. The Box-complex algorithm ' 


has been used for this 
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TABLE V 

STARTING ‘CONCENTRATIONS FOR SSP* 



INITIAL 

CONCENTRATIONS (mol/kg) FOR 

DIFFERENT 

RUNS OF SSP 



= 22 
n , o 

35 

55 

72 

95 

[ c i 3 


0.0000 

0 . 0000 

0.0000 

0.0000 

0.0000 

E p i3 


0.00075 

0 . 00068 

0 . 0006 

0.00055 

0.00049 

X 

o 


0.03122 

0.04331 

0.05244 

0.05456 

0.05367 

X 1 


0.69388 

1.53307 

2.93076 

3.95800 

5.13901 

A 2 


22.10 

80.43 

260.10 

485.20 

902.20 

[W] 


0.0000 

0.0000 

0.0000 

0.0000 

0.0000 


* Initial conditions in liquid phase polymerization to obtain the 
above values : 


Temperature 

= 

513 . 

0 K 

t C l'o 

s 

8.8 

mol/ kg 

^o 

as 

0.0 

mol/kg 

A 

o, o 

= 

0.0 

mol/kg 

X. 

1,0 

= 

0.0 

mol/kg 

X 2,o 

= 

0.0 

mol/kg 

[W] Q 

= 

0.16 

mol/kg 






purpose. The objective function that has been optimized in the 
program is 


PK^e^-i^L^^; D w ) = 


Max [- s abs(U n, e xptl»> ~ Vtheor<:»> } 
j=l ^n,theor^ 


where N is the total number of experimental data points for all 
M n 0 reported 5 , and the subscripts, exptl and theor, indicate the 
experimental and theoretical values of w n » 
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CHAPTER 3 

RESULTS AND DISCUSSION 

The computer program was tested for several special cases to 

ensure that it was free of errors. The number of finite difference 

grid points was increased from 11 to 21 and the results on M n (t) 

were found to be essentially the same (within 1% at t = 24 hr) . 

Similarly, the results were unaffected when the parameter, TOL, in 

the NAG library subroutine, D02EBF, was decreased from 10” 4 to 
—8 

10 . . After these tests, the Box complex method was used to obtain 

the best-fit values of the seven curve-fit parameters, using both 

Eqns. (17) and (18) for k. . Table VI gives the final values of 

1,0 

these parameters for both these cases. 

Simulation results for the five different experimental runs 5 

(corresponding to different values of the initial number average 

chain length, v n Q ) using the best-fit values of the parameters 

(Table VI) are shown by the solid curves in Figs. 2 and 3 for the 

case when the intrinsic rate constants, k, , are proportional to 

1,0 

the zeroth moment (Eqn. 17) . The experimental data points are 
also shown in these diagrams. Excellent agreement between theory 
and experimental data is observed. The dotted curves in these 
diagrams show the theoretical results using the best-fit values of 
the parameters with k^ Q - (Eqn. 18) . The final values of the 
objective function, F (in Eqn. 21), for these two choices of k. 

1 , U 

are -2.53 and -1.89 respectively. Based on these values as well as 
on the agreement between the theoretical curves and experimental 
data in Figs. 2 and 3 it is difficult to say which of these two 
correlations for the intrinsic rate constants is superior, and we 
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TABLE VI 

BEST-FIT VALUES. OF THE PARAMETERS USED FOR THE 
SIMULATION OF SSP OF NYLON 6 AT 190 °C 


Optimal Parameters 


k i,o = Cj 


k i,o = 


©! (hr) 
© 2 (hr) 

© 3 (hr) 


4.626x10 

406.557 


6476.71 


3.339X10 

486.251 


2586.66 


49.183 (kg 2 /mol 2 -hr) 0 . 3247 (kg/mol-hr) 
14996.9 (kg 2 /mol 2 -hr ) 405.7 (kg/mol-hr) 
360.59 (kg 2 /mol 2 -hr) 7 .7116 (kg/mol-hr) 


ID w (m /hr) : 

Other Parameters / Conditions 


= 2*10~ 4 m 


[W] = 0.00001 mol/kg 


3.366x10 


1.065x10 


= 463 K 
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have arbitrarily selected the form in Eqn. (17) for further study 
in this work. A comparison of the present results with those from 
our earlier work 10 which used the Fuj ita-Doolittle theory for the 
micro-level diffusivity reveals that the agreement of the present 
theory is better. In addition, the present theory is more elegant 
since the parameters of the Vrentas-Duda theory are obtained by 
independent measurements on nonreacting systems. 

The sensitivity of the results to changes in the parameters 

is now studied (using Eqn. 17 for k. ) . The value of u of 55 

I/O n / o 

is selected for this purpose. Figs. 4-8 show the effects of 
varying the seven different parameters one by one, about their 
best-fit values, keeping the other parameters at their reference 
(Table VI) values. The effect of 9 ^ on u n is shown in Fig. 4. 
Decreasing speeds up the reaction, particularly in the initial 
stages, and increases the final values of h n . Lower values of 
signify higher values of the micro-level diffusivity of the 
polymer molecule in the reaction mass. This leads to an increase 
in the forward rate constant of the polycondensation reaction 
(other quantities remaining unchanged) and hence in the values of 
U n . Very high values of 9 1 ( >0.01 hr) suppress the forward step 
of the polycondensation reaction so much that depolymerization 
takes place during the early stages. It is interesting to note 
that this figure illustrates how the presence of diffusional 
limitations (non zero values of 0^ ) reduces the value of 
significantly. 

The effect of varying the parameter, © 2 , on the progress of 
reaction is shown in Fig. 5. Increased values of 9^ lead to a 
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decrease in the micro-level diffusivity of the water molecule and 
a corresponding decrease in the reverse rate constant, (other 

quantities remaining unchanged) . This means that the 

polycondensation reaction is forced towards the right, leading to 
higher values of JI n . The sensitivity of the results to © 2 is less 
than that towards 9^. 

The results are found to be relatively insensitive to the 
value of © 3 . Decreasing © 3 from the reference values of 6476.7 hr 
to 100 hr decreases the value of jl at t = 24 hr only from about 
246.1 to 220.65. A similar insensitivity of is observed to the 
value of (between 10 hr to 1000 hr) . In fact, use of 0 

(which makes k. = 0) also gave the same u history. This means 

X f O XI 

that the ring opening reaction in Table I can easily be omitted. 
Indeed this reaction was omitted in our previous study} 0 based on 
intuitive reasoning. It may be emphasized that we cannot omit the 
polyaddition reaction from the kinetic scheme even though we get 
the same results with © 3 > » (i.e., k 3 = 0). This is because 

f 

the rate constant, k_ of the reverse step does not become zero 

J f o 

as © 3 » oo . Indeed, this reverse reaction is responsible for 

increasing the value of [C^] in the polymer particle from its 
initial value of zero, as the polymerization progresses. 

Fig. 6 shows the effect of varying C 2 * Increasing the value 

2 —1 —i 

of £ 2 beyond the reference value of 14996.9 to 22000. kg mol hr 

does not affect the progress of the reaction much, but lowering 

the values of does lead to considerable effect on the progress 

of the reaction. In fact, it is observed that as C 2 is increased 

2 —1 “1 

from a low value of 100 to 22000 kg mol hr , h n (t) goes up with 
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time more rapidly near t « 0 (due to higher rates of the forward 
step in the polycondensation reaction) , but flattens out at larger 
values of t to give lower final values of h n . The latter is 
because diffusional limitations become more pronounced and occur 
earlier for higher values of C 2 * Fig. 6, thus, illustrates the 
effect of competing physical phenomena. A similar phenomenon is 
observed in Fig. 7, where is varied (note C 3 going from 500 to 
1500 kg 2 mol“ 2 hr _1 ) . The curve for C 3 =0 (which makes k 3 = 0 also) 

in Fig. 7 shows what happens when we omit the polyaddition 
reaction from our kinetic scheme. This behavior is to be 
contrasted with the curve for £ 2 = 100 kg 2 mol 1 in Fig. 6. The 

importance of the polycondensation reaction over the polyaddition 
reaction is quite evident. 

The macro-level diffusivity of water through the reaction 

mass, E> w , (or the degree of crystallinity) plays an important 

role in SSP as shown in Fig. 8. For values of D w above about 

10 ~ 10 m 2 /hr (corresponding to higher fractions of amorphous 

material) the values of h n are higher. Values of ID W below the 

reference value of 3.37X10 -11 m 2 /hr lead to lower The results 

-13 

become relatively insensitive to below a value of about 10 

m^/hr. For such cases, the water produced by chemical reaction 

does not diffuse out much, and the reaction is driven towards 

final values of h n controlled by equilibrium alone. One need not 

solve partial differential equations for such cases. It should be 

emphasized that the degree of crystallinity would affect not only 

0 but also the values of 6. and its influence would be more than 
w A 

depicted by Fig. 8. 
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The effect of varying some important operating conditions has 

also been studied. The increase in ZZ is found to be relatively 

n 

insensitive to the surface concentration of water, [W] , when its 

s 

value is below about 10 ^ mol/kg (see Fig. 9) . For values of 

—3 

[W] s above about 10 mol/kg, however, the results are found to be 
sensitive to its value; This implies, technologically, that one 
needs to reduce the concentration of water in the vapor phase 
surrounding the nylon 6 to a threshold value, and not much is 
achieved in going below this level. It is interesting to note 
that the difference in [W] g shows up only after about 6 hr of SSP. 
The radius, R, of the nylon 6 powder (assumed spherical) is within 
the range of 0.2 mm - 0.5 mm, as reported by Gaymans et al . The 
reference value used in this study is 0.2 mm. It is found that 
(see Fig. 10) decrease in the value of R from 0.5 to 0.2 mm causes 
an increase in the final value of of about 40 to 50. It is 
interesting to note that much more is achieved by reducing the 
pellet size from 0.3 mm to 0.2 mm than from 0.5 mm to 0.4 mm. 
This is becuase the major part of polymerization occurs in a thin 
shell near the periphery (see later) . Figs. 9 and 10 quantify the 
effect of two of the most important operating conditions. 

We also carried out a simulation run for the case when the 
monomer and water are not leached out of the polymer particles 
before the SSP, using the reference values of all the parameters 
(for the case when k A Q - C* * c ) • The values of [C 1 ] Q and (W] Q are 
5.869 and 0.09683 mol/kg respectively (for M n , 0 = 55). A sudden 
jump in M n is observed near t - 0, but the values for t larger 
than about 9 hr, are much lower (213.7 as compared to 246.1 at t = 
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24 hr) . The sharp initial rise in ji n is due to high initial 
concentration of the monomer which drives the polyaddition 
reaction to the right-hand side. This effect dominates over the 
effect of higher [W] values, which drives the polycondensation 
reaction to the left (depolymerization) . The lower final values of 
U n emphasizes the need for carrying out leaching before SSP. 

Fig. 11 shows the results for SSP of nylon 6 with a single 

intermediate remelting (homogenizing) of the polymer particles. 

Reference values (Table VI) of the parameters have been used. 

Points p, q, r and s in this figure show that if the total time 

for SSP (i.e., before and after remelting) is the same (24 hr), 

higher 4 n product can be obtained by intermediate and early 

homogenization. The results are not as dramatic if the 

homogenization occurs late. This trend is qualitatively similar 

5 

to the experimental observations of Gaymans et al . The higher 
values of n n result from lower (uniform) values of [W] in the 
inside of the polymer particles, resulting from homogenization. 

Fig. 12 shows the water concentration profiles at different 
values of t (in absence of remelting)' for the reference 
conditions. It is observed that the water concentration builds 
up in the interior of the polymer particle (due to generation) 
quite early during the reaction. However, macro-level diffusion 
of W leads to a relatively sharp decrease in [W] in the outer 
shell of the particle, 0.16 mm < r < 0.2 mm. The lower water 
concentrations in the outer shell are associated with 
significantly higher values of #1^ in this region, as shown in Fig. 
13. This also explains why remelting of the polymer particle leads 
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intermediate homogenization (remelting) of polymer 
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55 and with Eqn. 17) . 




FIG. 13 

Radial profile of ji at different times for the 
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to higher product, as shown in Fig. 11, since intermediate 

homogenization removes the spatial gradient, and reduces the 
averaged value of [W] in the inner regions. It may be added that 
these profiles have been generated using 21 finite— difference 
grid points because of the relatively sharp gradients near r = R. 

Fig. 14 shows how the monomer concentration increases from 
its initial value of zero as the reaction progresses. Since it 
has already been found that the ring-opening reaction is 
unimportant, this means that the reverse step of the polyaddition 
reaction is responsible for the generation of the monomer. The 
reason why [C 1 ] is lower near the outer fringes of the polymer 

particle is because the total concentration of polymer (i.e., £ 

< 

[P n ]) molecules is lower there (associated with the higher values 
of 4 n ) . 

The variation of with location would be expected to lead 
to values of the spatially averaged (mixed) polydispersity index, 
Q, above the normal value of 2.0. Fig. 15 shows how Q increases 
with time for the five cases studied. Q is found to rise 
rapidly at first, and then increases more gradually to values 
above 2.0. The initial values of Q are not 2.0 for the low a 

n, o 

cases due to low conversions at the beginning. 

It is interesting to note (see Fig. 16) that the apparent 
rate constant, k 2 , increases with time, rather than decreases. 
This diagram shows k 2 at two positions, one near the center of the 
spherical polymer particle (at r = 0.04 mm), and the other near 
the periphery (at r = 0.18 mm) for the reference run (with u = 
55) . The value of the intrinsic rate constant, k 2 o' found to 
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FIG. 14 

profile of [C ] at different times for the 
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decrease with time at both these locations. The latter is 

associated with a drop in the value of X Q (see Eq. 17) which, in 

turn, is due to the increase in the value of u n with time. Fig. 13 

shows that 4 is higher at r = 0.18 mm, and so X and k_ would 

0 2/0 

be expected to be lower at this location than at r = 0.04 mm. This 

is indeed confirmed in Fig. 16. Fig. 17 shows that the ratio, 

^2^2, o' i ncr ©ases w ith time at both locations. This increase 

indicates that the diffusional resistance decreases as the SSP 

progresses. This behavior is in sharp contrast to the increase in 

the diffusional resistance in the polymerization of methyl 

methacrylate (gel effect) . The decreasing diffusional resistance 

in the SSP of nylon 6 is associated with the decrease in the 

volume fraction, 0 , of the polymer with time due to generation of 

C 1 and W. In fact, Figs. 12, 14 and 17 show that as t increases, 

[C^] increases at both locations but [W] increases continuously at 

r = 0.04 mm only ([W] decreases with time at r = 0.18 mm after an 

intial rise) . This leads to <p^ being higher at r = 0.18 mm. This 

complex interplay of several phenomena — generation of and W, 

the macro-diffusion of W, decrease of X and with 

o 2,0 

polymerization, etc., in influencing the rate constant, k 2 , is 
clearly brought out in these figures. Similar increases in k_/k 

J O / u 

and k'/ki with time are also observed. 

2 ' 2,o 

It may be mentioned here that most of the simulation results 
reported herein have been generated using the closure conditions 
given in Table IV. Use of [P 2 ] = [P-J gives negative monomer 

concentrations for some extreme choices of parameters which are 
far removed from the reference values, e.g., for very low and 
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very high Whenever this is encountered, we have used the 
closure condition, = 0*001 [P^] . The moment closure equation 
in Table IV is used unchanged. Results were regenerated using this 
new closure equation for all the other parameter values too and 
were found to be insensitive to this change. We recommend the use 
of this new closure condition in future studies. 



CHAPTER 4 

CONCLUSIONS 
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A new model having a strong molecular basis has been proposed to 
account for the effects of segmental diffusion on reversible step 
growth polymerizations. This model is more fundamental and general 
than earlier models on SSP and can be used to analyze a variety of 
SSP processes, as for example, nylon 6 polyethylene terephthalate, 
polybutylene terephthalate, nylon 66, etc. The model is applied to 
the specific case of the SSP of nylon 6 in this study. A comparison 
of the present results with our earlier work 10 reveals that agreement 
of this theory with available experimental data at 190°C is far 
better. 

The optimal values of six parameters, and ; i = 1,2,3 (or 

five, with = 0 and © 3 > «) obtained in this study can be used at 

190°C (alongwith the model equations and the new closure condition 
developed in this work) to design and optimize industrial nylon 6 SSP 
reactors in the future. The seventh parameter, ID w , used in this 
study depends on the crystallinity of the sample, and must be 
estimated using pilot-plant studies, till such time that a better 
understanding relating it to the morphology (which could be 
independently studied) becomes available. It may be added that 
similar optimal-parameter estimation could be performed for SSP of 
other polymers (or for nylon 6 sit other temperatures) , to give a more 
rational design methodology for those reactors as well . This study 
quantifies the effects of three important operation parameters, 
namely, radius of the polymer particle, water concentration in the 
vapor phase, and the degree of crystallinity of nylon 6 on the 
spatially-averaged molecular weight history. 
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INPUT FOR OPTIMIZATION OF PARAMETERS 

I. BOX-COMPLEX ALGORITHM PARAMETERS 

a = 1.3 

|3 = 0.0001 
7 = 5 

5 = 0.00001 
IC - 1 
IPRINT = 1 
N = 7 
M = 7 
K = 8 


II. FOR THE FORM k. ^ = r. A OF INTRINSIC RATECONSTANTS 

1,0 1 o 

1* Initial Guess Values of Parameters 
0^ = 3.25 X 10" 4 
e 2 = 380.0 
0 3 = 1475.0 

= 7.25 X 10" 2 
C 2 = 4150.0 
C 3 = 135.0 

D r = 4.3 x 10 -11 m 2 /hr 
w 


2. Range of Parameters 


3 

<1 

^2 

£3 

D w (m 2 /hr) 


Lower bound 
l.Oxlo' 4 

5.0 

10.0 

1.0xl0~ 2 

l.OxlO 3 

5.0 

11 

1.0x10 


Upper bound 
6.0X10 -4 

1000.0 

10000.0 

100.0 
20X10 3 
800.0 

11 

5.0x10 


3. Final Error (%) 


3.78 



III. FOR THE FORM k. = £. OF INTRINSIC RATECONSTANTS 

1.0 l 


1. Initial Guess Values of Parameters 

0 1 = 4.65 X io " 4 

0 2 = 250.0 

0 3 = 635.0 
= 0.0113 

C 2 = 331.0 
C 3 ~ 13.5 

D w - 3.4 x 10 -11 m 2 /hr 


2. Range of Parameters 

0 1 

0 2 

0 3 

^1 

<2 

£3 

D w (m 2 /hr) 


Lower bound 
l.OxlO -4 

5.0 

10.0 
0.001 
10 
1.0 

l.OxlO 11 


3. Final Error (%) = 5.06 


Upper bound 
6. OxlO -4 

1000.0 

5000.0 

10.0 

1000.0 

50.0 

11 

5. 0X10 J ‘ X 


CENTRAL 

4* jfo A 
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0F SOL!D STATE POLYMERIZATION* 


-CMUKUND R. KULKARNI 92102123 

,« B .« a « a « 8 „ #ssss# ; BSBaeaas 

Using dOEebfe.f routine fro* 


OF NYLON 6 


ssiasaKSSSE; 

few 


NAG 1 i brarv 

sattwajsaasseaaaaaaaa^aias'-*-.-.^- ** 

ASPECTS : ==«=»==*=:=:«==,= 

ALL THREE REACTIONS CONSIDERED 

HACRO-LEVEL DIFFUS.VITV OF UATER IS A PARAMETER . 
0.001 CPU is THE CLOSURE CONDITION. 

1 for k i , o * zetai*Ia*o 

2 for k i , o = zet a i 
points np a ti. 


CP23 

USE icase a 
lease a 
5. No. of grid 
INPUT FILES ; 

' input ' Give 


initial conditions 


9ive icase 4 values of the parameters 

IMPLICIT DOUBLE PRECISION (A-H,Q-Z) 

DOUBLE PRECISION N INDEX 
. . Parameters . , 


PARAMETER 
PARAMETER 
PARAMETER 
. . Local Arrays 
DOUBLE PRECISION Udll), 

.. External Functions 
EXTERNAL D02EJW, G 

.. External Subroutines .. 
EXTERNAL D02EJY, 

EXTERNAL D02EBF 

. . intrinsic Functions . 
INTRINSIC DBLE 


(N«66, IW=21*N+23,NP=1 1 ) 
< N=6& , IWaSS00,NP*1 1 I 
(HOUTafi) 


Y(N) 


FCN, OUT, PEDERV 


DIMENSION STATEMENTS 

DIMENSION AO<3>,EOt3),AC(3>,EC(3),DH<3>,DSC3> 

COMMON STATEMENTS 
COMMON XEND, H, I 
COMMON /B1/Y 

COMMON /82/R CNpi, FKS2 (Np>, FKS3 C Np ) , RKS2 (Np), RKS3 C Np ) 

COMMON /B3/EK8, EK3 
COMMON /B4/ZETA2 , ZETA3 
COMMON /BS/DU, TEMP, MS, DR 
COMMON /BS/RR,PI 

COMMON /B9/THETA , AT , BT , TG INF IN I TY , P 1 , ALPHAL , ALPHAS 

COMMON /B1 1 /RADIUS 

COMMON /B1 3/TGC NP ) > VFINP I , DD0 ( NP } 

COMMON /B 14/FX ! C NP ) ,FX8 C NP ) # FX3CNP > *FX4< Np ) , FXSC Np J , FXS C NP) 
COMMON /B 1 5/AVEC 1 , AVEP1 , AVELAMO , A VELAM 1 , AVELAM2 , UVEI# 

COMMON /B 1 6/ AWfe VP *, AVBT6 :• '.'V;.'' ' 

COMMON /B18/VMS,¥PS»SAmMa,VM0,2ETA,VFM ^7 

COMMON /BlS/NINDiX*MMrMEOt w / 

'COMMON /BSO/THETAHTMEf AS, YHETA3 7. 7 - "V. ,y , \-n ,• 'V- 



COMMON /821/FKS1 (NP),RKS1 <NP) 

COMMON /Bae/EKI , ZETA 1 
COMMON /B23/DD02 ( np ) , 0DO3 ( np ) 

COMMON /B24/VS0 , VFS, V88 
COMMON /B25/ZETA13, 2ETA23 
common /const/const 
common /case/icase 

OPEN STATEMENTS 

OPENCUNIT « 12, FILE * 'Taietal') 

OPEN ( UN I T » 14, FILE - 'input') 

OPENCUNIT * 15, FILE * 'input!') 

OPEN ( UN I T - 20, FILE « 'output') 

OPEN (UNIT - 21, FILE = 'ssp.oi') 

OPEN ( UN I T * 22, FILE - 'ssp.o2' ) 

OPEN (UNIT = 23, FILE * 'pdi') 

OPENCUNIT * 84, FILE * 'rateconst'I 

Executable Statements .. 

READING INPUT DATA 

DO I J «* 1,3 

READC 12,*)AO(IJ),EO(IJ),AC(IJ),EC(IJ), DHC I J ) , DSC I J ) 
END DO 

READ C 1 4 , * ) Y1 O,Y2O,Y30,Y4O,Y5O, Y60 

CONSTANTS 

RR : GAS CONSTANT J/K/gmoI 

Nav : AVA. NO. 

PI : 3.1417 
RR » 8 313D0 

....BOUNDARY CONDITION 

WS * 6.0000100 

INITIAL CONDITIONS .... 

DO JJ » 1 ,Np 
YiJJ) « YtO 
YCNp+JJ> * Y20 

Y < 2*NP + J J ) * Y30 
YI3*NP*JJ) * Y40 
Y(4*NP+JJ) « YSO 

Y C 5*NP+<J J ) « YSO 
IFC JJf .EO.NPITHEN 
YtS*NP+JJ) = WS 
END IF 
END DO 

TEMP » 463 1600 
TGINFINITY * 350.606 
Pt *"■ 6,0306 ■ v : ... . 

alphal *1.60-3 

ALPHAS :.4 v50-4-:v 4 vV •? ' , : - : ' 

RADIUS 0 OEtir 
NIN0EK«0;0<*0 



const * 0 . 00 1 dO 

reading parameter VALUES 

read! 1 5 , * ) i case 

read! 15,*) THETA 1 , THETAE, THETA3, ZETA1 , ZETA2, ZETA3, DW, TOL 
writeCEl , * ) 'THETA 1 , THETA2, THETA3,ZETA1 , ZETA2, ZETA3, DW, TOL ' 
writeiei , * ) THETA 1 , THET A2 , THEY A3 , ZETA 1 , ZETAE , ZETA3 , DW , TOL 
WRITE (20,331 ) 

1 format (lx, 'TIME ' , fix , ' AVEC 1 ' , fix , ' AVEP 1 ' , 4x, 1 AVELAMO' , 3x, 

+ ' AVELAH 1 ' ,£x, * AVELAM2 ' 

+ 4x , ' AVEW ' , 

+ fix , ' AVEMUN * ) 

FREE VOLUME PARAMETERS- 

VMS * 0 .86426300 
VPS * 0.81 232D0 
VS3 = 0 . 96909 
ZETA * 0 . 709£96d0 
GAMMA - 1 . 0000 
VMO *1.1 3E47D© 

VFM » 0.31 03300 
VSO * 1 .142900 
VFS * 0. 1795d0 

UM * 1 13. 16D0 
WSOL * 18.000 
ZETA13 * VMS/VPS* ( E . 0/3 . 0 ) 

ZETAS3 * ( VSS/VPS ) * ( W8/WH ) * ( 2 . 0/3 . 0 ) 

ZETA13 * 0.709E96D0 
ZETAE3 * 0 . 126500 

X * 0 . 0D0 

EQUILIBRIUM CONSTANTS 

EK 1 ■ 2.3516980-03 
EKE * DEXPC (0S(E)-DH(2)/TEHP)/RR) 

EK3 * D£XP( (DSC 3)-DH(3)/TEMP)/RR> 

EKE * t 033 . 24?4d0 
EK3 * £ . 456d0 


DR * O.OOEdO 


DO JK * 1 , NP 


R( JK) * DR* C JK - t ) 


END DO 



WRITE CNOUT,*) NYLON 6 SSP SIMULATION ' 

XENO *24.000 
MPEO * . 6 . ■ . 

ir * s V\ "■ '• ■ ’ ' 

WRITE CNOUT, ftstsi ' €alcal4ti®« with TOL TOL 

h ' *■ (XgN&^X i/D»LEf "■< I -;,-/";' ■ " 

' IFAIL * 0 . ■■■- • 

GALL BOSEBFlXVXENOjN, Yjt&L, IR,FCN,HPEB,jB02EJY,OOT, W, IW, IFAILT 
IF t TOL . LT. 0 .0100) WRI TE ( NOUT ,*) ■*. Range % ©o short for TOL * 



STOP 


99999 

99998 

99997 

99996 


FORMAT 

FORMAT 

FORMAT 

FORMAT 

END 


(IX, A, 08. 1 ) 

( 1 X , A , F7 . 3 ) 
(1X,A,6((1tE16.5)/)) 
( F4 . 2/ 6 ( ( 1 1E16.S)/) ) 


SUBROUTINE FCN ( T , Y , F ) 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

DOUBLE PRECISION NINDEX 
PARAMETER <N=»66,Np-1! ) 

DIMENSION F(N), Y(N) 

~ X / M N ?EH^ISr ,FKS3tNP1 ' RKSStNP, ' RK83<N ^ 

COMMON /BEI/FKS1 ( NP ) , RKS1 CNP) 
cottMon /const/const 
, . Executable Statements . . 

DR = 0 . 002D0 


DO £5 JJ » !,Np 

Uf * JJ 

ue « 1 »Np + JJ 

IJ3 a 8*Np + JJ 

IJ4 a 3*Np + JJ 

IJS = 4*Np + JJ 

IJ6 « 5*Np + JJ 

CALL RATECQNST (JJ,IJ3,IJ4, IJ6, Y) 


ZZZ a 
XL3 


F( IJ1 1 » 


F< IJ2) = 


F( I J3 ) 


FU J4> 


F( I JS 


IFCJJ 
YP . «*, 


* +FKSI (JJ)*YCIJt ) *YC IJ6) — RKS1 < j j ) *Y( I J£ ) 

“ Y(IJ5>*< 2 . 6D8*Y { I J5)*Y ( I J3) - Y( I J4 J*YC I J4) 
/ l Y ( I J4 ) *YC I J3 ) ) 

1 83 -FKSjt JJ>*Y( IJI J*YC IJ3J 

+RK83( JJ)*(Y(IJ3) - YIIJSJJ - ZZZ 

1 ° ~2 . 00O*FKSE( JJ)*Y( I J£)*Y< I J3) 

+£. 0D0*RK8£UJ}*Y< IJ6>*CYt I J3J - Y(IJ2) ) 

-FK83 (JJ)*YCIJ£3*Y{ IJI J 
+RKS3( J J ) *const *Y 1 IJ£) + ZZZ 

* -FKS£C JJT*Y(IJ3)*¥( IJ3) 

+RKSEC JJ)*Y( IJ6)*(YCIJ4) - YC I J31 > 

♦ZZZ ■ . 

« FKS3C JJJ*Y( IJt >*YC IJ3> - 
RKS3( JJ)*(Y(IJ3) - Y ( I J£ ) > 

+ZZZ . 

i> * g; 0D0*FKS2( JJ)4Y(IJ4)*YUJ4J 

+1 .9OS/3.0Da*RKS£(JJ)*Y(IJ6)*(Y(IJ4) - XL3> 

J 4^YtIJl >*«Yi I J3I + £ . 0D0*Y( I J4> ) 
t'Mti I J31 - £ . 9DD*Y( I J4 ) + YC IJ2) } 

■ +zzz :; iw> /;;• ■" ; 

.EQ.lj THEN 

.Y CRT ■■■■■ ' V ; ,V . -V ; , - , 



FCIJ6) * -Ft I J3> 

* +Dw*( (Y( IJ6+1 ) - 2 . ODO*Y( IJ6)+YP>/(DR*DRJ ) 
ELSE 

YP * Y ( i J 6~ 1 ) 

F< I J6 ) * -F ( I J3 ) 

* +Dw*C ( Y( IJ6 + 1 )-e. 0DO*Y( I J6 ) +YP ) /( DR*DR ) ) 

* +Dw*( Y( IJ6 + 1 } - YP)/R( JJJ/DR 
endi f 

F ( 66 ) * 0. OdO 
g5 CONTINUE 


RETURN 

END 

g 5 {BB«»a#a?aaaas«fl!»aasa!aaftisssa!saiaaa»!aaasassasasaasasBaaaaaBfc!assaaBSKasa 
SUBROUTINE 0UTCX,Y> 

IMPLICIT DOUBLE PRECISION CA-H,0-Z> 

DOUBLE PRECISION NINDEX 

* . , Parameters . . 

PARAMETER <N«66,NP-1 1 ) 

PARAMETER (NOUT*6) 

* . . Array Arguments . . 

DOUBLE PRECISION Y C N ) 

DOUBLE PRECISION RRR ( NP i 

COMMON XEND, H, I 

COMMON /B£/R ( Np ) , FKSS ( Np ) , FKS3( Np i , RKS£ ( Mp J , RKS3C Np ) 

COMMON /Bit /RADIUS 

COMMON /B 1 3/TC { NP ) , VF (NP t , DDO { NP > , AMUNCNP I 

COMMON /B14/FX1 CNF), FXE ( NP ) , FX3 C NP > , FX4 C NP ) , FX5 C NP > , FX6 C NP ) 
COMMON /B1S/AVEC1 , AVEP1 , AVELAMQ, AVELAM1 , AVELAH2, AVEW 
COMMON /Bt 6/AVEVF, AVETG 
COMMON /B17/FX7CNP) ,FX8(NPJ,FX9CNP) 

COMMON /B21/FKS1 CNP),RKSHNP) 

COMMON /B4/ZETAS, ZETA3 

COMMON /B20/THETAI , THETAS, THETA3 

UR I TEC 22, * > 'TIME®* , X 


do jjt 

» 1 

, NP 



Ui 

as 

JJ 



ue 

as 

t*Np 

♦ 

JJ 

IJ3 

as 

2*Np 

' 

JJ 

IJ4 

as 

3*Np 

* 

JJ 

IJ5 

as 

4*Np 

‘ + 

JJ 

US 

. as 

S*Np 

t 

JJ 


RRRCJJ) * RCJJI*RCJJI 
CALL R ATECOfiSf C J J , IJ3, IJ4, IJ6,Y> 
FX1 t JJJ a mnt Jj J*Tt Ut ) 

FX2C JJ )'-•* RRRIilJ J*YC I32J 
PX3 um ■ R R R tj J> * Y t I J 3 ) 
FX4tJjri * RRIHJJJ*Yf U4J . 
FXStJJI - RiR( W*Y( IJ5) , 
FX6 I JJ) « RRtUJJ)*YCIJ6 ) 

■: FX7CJJ) « AHUNCJJ) ■ ■ ■ 



i'i FX8 < J J ) » RRR ( JJ ) *VF( JJ ) 

FX9 1 J J ) * RRR( JJ)*T6( JJ) 

WRITEI20, 19>X,Y(IJ1 ),Y(IJ2),Y<IJ3),Y(IJ4> 

+ , Y( IJ5) , Yt IJ6) , AMUNt JJ) 

19 FORMAT (8(E9.3 f 1X)) 

URiTEtee, *)R( JJ) , ' aaunC j j ) 

END DO 

CALL AVERAGE! Y, AVEMUN) 

avepdl » avela®E*avela»o/(avelaat**2) 

WRITE<6, *) 'TIME** ',X,' AVEMUN : AVEMUN PDI : ' ,av«pdi 

WR I TE ( 20 , £9 ) X , AVEC1 , AVEPt , AVELAMQ, AVELAM1 » AVELAM2 , A VEU , AVEMUN 
29 FORMAT (8(E9.3,1X)) 

WR I TE <21, 22 ) X , AVEMUN 
WRITE ( 23 , 22 ) X , avepdi 
22 FORNAT(4X,FtO.3»4X,F10.4) 

X * XEND - DBLEt I )*H 
1 = 1-1 
; RETURN 

* '' 

99999 FORMAT C F 1 9 . 4/6 ( f 6F1 6 . 8/5FI 6 . 8 ) / > ) 

END 


SUBROUTINE RATEC0M8T ( JJ , IJ3, IJ4, I J6 f Yl 
IMPLICIT DOUBLE PRECISION tA-H,0-Z) 

DOUBLE PRECISION NINDEX 
PARAMETER CN=66,Np=1 1 ) 

DIMENSION Y < N ) , FKSEOC NP ) , FKS30 C NP ) , RKS29 ( NP ) 

COMMON /B2/R (Np) , FKS2 CNp) , FKS3 C Np ) , RKS2 C Np ) ,RKS3 I Np > 
COMMON /B3/EK2, EK3 
COMMON /B4/ZETA2/ ZETA3 

COMMON /B9/THETA, AT,BT , T6INFINITY, PI , ALPHAL , ALPHAG 
COMMON /B1 3/TStNP ) , VFCNP ) , DDQINP ) 

COMMON /B20/TMETA1 , THETA2, THETA3 
COMMON /B21/FKS1 <NP),RK81 (NP) 

COMMON /B22/EK1 ,ZETA1 
COMMON /B23/DD02I np ) , DD03 C np ) 

COMMON /B25/ZET A f 3 # ZET A23 
common /case/icase 
GALL BIFFUSIVITYUY, JJ,i }3, i j4»i j6) 

XYZ * THETA/ BOOT JJ) 

XY21 * THETA1/DD0C JJ) 

XY22 * THET A2/BBBS f J J ) 

XYZ3 * THETA$/Bt®3CJJ) 


FKSUJJI * 2Emi 
FKS2,0tJJJ^*-ZltAe 
FKS30I JJ) « ZETA3 




i f ( i case . ne . 1 ) go to £££ 
FKS1 i JJ) » ZETAf *Y( IJ3) 
FKS20CJJ) * ZETA£*Y( IJ3) 
FK830 C J J ) - ZETA3*Y(IJ3> 


88fi RKSl(JJ) ■ FK8t ( JJ)/EK1 

RKS£0<JJ) - FK8E0 ( JJ ) /EKE 
RKS3 ( JJ ) « FKS3Q ( J J ) /EK3 

C diffusion affected rate constants 

FKSE(JJ) = 1 . 0D0/( 1 , 0D0/( FKSEO < JJ ) ) + XYZ1*Y(IJ3)) 

RKSE(JJ) = 1 . OD0/ ( 1 . 0DO/( RKSE0< J J) ) + XYZ8*Y(IJ3>) 

FKS3CJJ) » 1 .000/(1 . 000/ ( FKS30 < J J ) ) + XYZ3*Y(IJ3)) 

RETURN 

END 

'' jJsss5 5:ss:a3;ssi;ss:!ssss!E3;c5as5ss!ssas,'sasiasassc:saatsa3Ba#asBssBBfiBaESBSSB!#i!SSSBS8ES#i 

! SUBROUTINE D I FFUSI VITY1 < Y , j j , i j 3 , i j 4 , i j 6 > 

IMPLICIT DOUBLE PRECISION (A-H,0-Z> 

DOUBLE PRECISION NINDEX 
PARAMETER ( N*66 , NP*1 f > 

DIMENSION Y(N j ,DD(NP J ,DD1 (NP),DD£CNP) 

COMMON /BS/DM, TEMP , US , DR 

COMMON /B 9 /THETA, AT , BT , TC INFINITY , PI , ALPHAL, ALPHAS 
I COMMON /Bit /RADIUS 

COMMON /B 1 3/TG ( NP > , VF < NP > , DD0 ( NP ) , AHUN C HP > 

I: COMMON /B t 8 /VMS , VPS , GAMMA, VM0 , ZETA , VFM 

COMMON /B1 9/NINDEX, WM,USOL 
* COMMON /BE0/THETA1 

COMMON /BS3/DD0E ( np > , DD03 < np ) 

COMMON /BS4/VS0, VFS, VSS 
COMMON /BE5/ZETA1 3, ZETAE3 

AMUN( JJ) * Yf IJ4)/Y( I J3) 

TG ( J J ) - 1 . 000/ ( I . OD0/T6 INFINITY + Pf/AHUN(JJ)> 

VF ( J J ) * O.0ES + (ALPHAL - ALPHAS ) * ( TEMP - T6CJJJ) 

DD0CJJ) * DEXPtE. 303*VF( JJJ/CAT + BT*VF(JJ)I» 

VPO ■ VPS*(t . 0+VFC JJ) ) 

TV *» (YC JJ)*MW*VH0 + YC IJ3)*AMUN( JJ)*WM*VPO + YC IJ6 )*WSOL*VSO) 
PHIM » Y{ JJI*«M*VH0/ TV 
PHIS » Y ( I J6 ) *WSOL* VS0/TV 
PHlP » YUJ3i*AMUN( JJ)*WM*VP6/TV 
' OH * 1 . 0/VM0 . 

■ . DP « 1.0/VPO \ 

DS « 1 . O/V60 

wP*.-V - VFP »■. VF( J J f ' ' ' • • ; 

i XNUM ^ SA«MAW<0M*PHIM*VMS/ZETAf3 + D3*PHIS*V8S/ZETAS3 ♦ 

Iv-'.':.- * r ■ ..DP*PHI#.*VP^Ift> . • ; V .'^'p 

I; XDEN = DM»MIW*VM#*VFM * DS*PHIS*VSS*VFS *BP*PHlP4VPSeVFP 

chi-' «•' XNUtf/XDfEit 'Yi : \.#ov 1 j? 

^;-- v . 'ch x ^ '-''ft V v' ■ / v ■ ■ ■ £ ; 


DO ( JJ ) ■ DEXP ( -CHI + CHIO) 

DD02CJJ) » DEXP t ZETA23* t-CHI+CHIO) ) 

DD03( J J ) * DEXP CZETA1 3* ! -CHI+CHI 0 ) ) 

DD© < J J ) » !AMUN(JJ)**!~N INDEX) ) * DD(JJ) 

RETURN 

END 


SUBROUTINE AVERAGE ! Y , AVEMUN ) 

IMPLICIT DOUBLE PRECI SION ( A-H , O-Z ) 

DOUBLE PRECISION NINDEX 
PARAMETER !N®66,NP=1 1 ) 

DIMENSION YIN) 

COMMON /B11 /RADIUS 

COMMON /B13/TG(NP) ,VF!NP> ,DD0!NP> .AHUNtNP) 

COMMON /B1 4/FX1 (NP ) ,FX2(NP ) ,FX3!NP ) , FX4 C NP > , FX5( HP ) , FX6 1 NP ) 

COMMON /B17/FX7(NP) ,FX8(NP) ,FX9(NP ) 

COMMON /B15/AVEC1 , AVEP 1 , AVELAMO , AVELAM1 , A VEL AMS , AVEW 
COMMON /B1 6/AVEVF , AVETG 

DENO » 10.08**3.0 )/3 . 0D0 

CALL INTEGRATION! FXt ,0.0,0. 08, 1 1 , AINT1 ) 

CALL INTEGRATION! FXE ,0.0,9.08,1! , AINT2 ) 

CALL INTEGRATI0N!FX3, 0 .0,0.02,11 , AINT3) 

CALL INTEGRAT10N!FX4, 0 .0,0.02,11 , AINT4) 

CALL INTEGRATI0N!FX5, 0.0, 0.02,11 ,AINT5) 

CALL INTEGRATION! FX6, 0.0,0. 02, 1 1 , AINT6 ) 

CALL INTEGRATION! FXT ,0.0,0.02, 1 1 ,AINT7) 

CALL INTEGRATION! FX8 ,0.0,0.02, 1 1 , AINT8 ) 

CALL INTEGRATION! FX9 ,0.0,0.02,11, A I NTS ) 

AVEC1 » AINT1/DEN0 
AVEP 1 * AINT2/DEN0 

AVELAMO « AINT3/DEN0 
AVELAM1 « AINT4/DEN0 
AVELAM2 * AINT5/DEN0 
AVEW «= AINT6/DEN0 

AVEVF * AINT8/DENG 
AVETG » AINT9/DER0 
AVEMUN * AINT4/AINT3 

RETURN 

END 

«sas>ssss:8ssaxaeB3aBssss:saaeassissssaB3SB3sasssas«sksBiesssais*B«seBa 

SUBROUTINE INTEGRATION! FX, XI , XL, NPOINTS, XINT) 

THIS SUBROUTINE INTEGRATES FUNCTION 'FX* OVER THE RANGE *XI* TO *XL' 
' NP0INT8-1 * NO. OF INTERVALS. NP0INT8 SHOUD BE ODDI 
* XIMT * IS THE VALUE OF INTEGRAL USING SIMPSONS* RULE. 

IMPLICIT DOUBLE PRECISION <A-H,0-Z) 

PARAMETER /. * i 

DIMENSION FXINPJ ' 

: ■ mV *:NP0INTS , ■ 

\M£^' npojnts 



DX * O.O08 


sum *0.0 

DO JJ * 3, NS, £ 

sum * sum + fx( jj> 

END DO 

SUMS * 0.0 
DO JJ * £,N1,S 
SUMS * SUMS + FXCJJJ 
END DO 

SUM * FX( t ) + FXC NPOINTS ) + S.0D0*SUMi + 4. 0D0*8UM£ 

XINT * DX*SUH/3 . O 

RETURN 

END 

^saa*aa*#»»J5EB#s##saaisasss«ssassB:asasassas : ssssssssBa«»«8* ! =* s 



OPTIMIZATION USING BOX-COMPLEX ALGORITHM 


************ ***********m*******^t****m***m******^*****m*********** 
MAIN LINE PROGRAM FOR COMPLEX ALGORITHM OF BOX 
*************************** 

IMPLICIT DOUBLE PRECISION A-H, O-Z ) 

PARAMETER ( N*7, M«7 , K«8 > 

Double precision X( K, M J , FC K } , G(M> , H( M > , XC( N> 

INTEGER GAMMA 
COMMON /B27/YINI (5,6) 

COMMON /BE8/EXMUN(5, 13) ,NSET 
COMMON /B3Q/ERR0R, TOL 
COMMON /B31/IT, ITHAX 
common /bb/ichoice 
common /const/const 
OPEN STATEMENTS 
OPEN! UNIT « 14, FILE » 

OPEN( UNIT = 33, FILE « 

OPENCUNIT * El , FILE * 

OPEN ( UN I T*S5» FILE*' guess ' > 

OPEN ( UNI T®59 , FILE*' opt * } 

OPEN(UNIT«60,FILE*'QtttO©?' > 

N0=»59 

INPUT FILES : 

'input' : Initial Concentrations for Ail the Sets 
'exdata' : Experimental Data (5 Sets! 

'opt' : output 

t«4*t**t«tti|i*4tiktt«t*t**t***«*f*******t********t***************t*4 
const “ 0 , 001 dQ 

****** * ******^ ****^*** *^**9 ******** » ****** *****mrn ********** **4 


' input ' ) 
'exdata' 1 
ssp.ot ' ) 


DATA IC, IPRINT/f ,1/ 

BATA ALPHA, BETA, GAMMA, DELTA/1 .3,0.0001,5,0, ©00001/ 
read( 55, * ) i choice 
if C ichoice .eq. ©)go to 1EE 
do i j k * 1 , K 

READ (55, *)<X{ijk,JJ),JJ=1,M) 
end do 
go to 558 

2 READ( 55>*HX(t,JJ),JJ=1,N) 

8 READC 14, *)nset,TOL» 1 TMAX 

do 77 iset =1 , nset 
read! 14, * M YINIC is#* , J k J , jfc*t ,6) 
read( 33, ♦ ) (EXKUNt iset ,ij),ij*1,13) 

■ continue . 

************************************* ************n 

■ WRITE(NO,010> " . ' 

F0tMAT<//, i8X, 'COflPtEX PROCEDURE QF BOX' ) 

FORMA 

' WRITEiNO* ®t t lit, II, Xii tMAX# 1C, ALPHA, !E'fA;*6AMHA, DELTA 
FORMATC //jEX, 'N® 'jl4, 3X, 'H*', 14, 3X, 'K*' ,14, EX, MM AX*' , 



1 

a 


16, ex, ' IO' , l4 ( //,ex, * ALPHA- ' , 
' GAMMA 8 ' , 14, 3X, ' DELTA® ' , Ft 0 . 5 ) 


F8 . 4 , 5X , 


BETA® 


Ft 0 .5,3X, 


_ CONSX(N,H,K, ALPttA, BETA, 6AHHA, DELTA, X,F, IEV2, 

3 NO i6iH,XC, IPRIHT ) 

IF ( I T- 1 THAX >20,20,30 
MR1 TE ( NO , 0 t 4 )F( IEV£) 

FORMAT! ///, 6X, 'FINAL VALUES OF THE FUNCTION ®',E20.8J 
MR I TE C NO , 0 1 S ) 

FORMAT { /// , ex, 'FINAL X VALUES') 

DO 300 J=1 ,N 

WRITE! NO , * ) J.XCIEVE, J) 

6 FORMAT!/ ,8X, ' J*' , 14, SX, *X=' , EE0.8) 

CONTINUE 
GO TO 909 

MR I TE ( NO ,017) ITMAX 

FORMAT!///, EX, 'THE NUMBER OF ITERATIONS HAS EXCEEDED', 14 
6 'PROGRAM TERMINATED') 

STOP 

END 

** + *#**m*#*##mmm#+**#* **********m*******m***m**m************* 
SUBROUTINE CHECK t N, M, K , X , G , H, I , KGDE, XC , DELTA, K I ) 

IMPLICIT DOUBLE PRECISION! A-H, 0-Z> 

Double precision XtK,M) ,G!M) ,H!M) ,XCCN) 

INTEGER GAMMA 


KT*0 

CALL CONST! N, M, K, X, G, H, I ) 


CHECK AGAINST EXPLICIT CONSTRAINTS 


DO SO J= 1 , N 

IFCXU , J)-G! J) )E0,e0,30 
X! I , J )«=G ( J J+DELTA 
GO TO SO 

IF!H( J)-Xt I, J))40,40,50 
XC I , JT )-H! J ) -DELTA 
CONTINUE 

IFCKGDEJt 1 0,t 16,60 
NN«N+t 

DO 100 J *NN, M 

CALL CONST! N,H,K,X,G,H, I > 

IF!X( I, J)-CCJM80,70,70 
IF ( H ( J J -X! I, J>) 86, 100, tOO 
lEVt “I . 

KT*t ■ 

CALL CENTS (N,M,K, IEVt , I ,XC,X»K1 J 
DO 99 JJ »t ,H 

X<I, X C I , JJT+XC l JFJ 1 1 /SO 

continue ; : 
coNTiiMiE 'v\7\h7:7:./i'7:v .■ •" ■ • . 
IFLK? ntov : :7^77 
RETURN ■. -:77V: > ■>'- V-r ^"7 ' ; - 

•: :- o:.-.7^ ..,7 ; .7-7:7: ; ' 7:.... ' 7v'7 



' ' 1 > ***** lf *****^ ;,! **^***%*^*#********^****;***##***##************ 

SUBROUTINE CENTR(N,H,K, IEV1 . I ,XC,X,K1 > 

IMPLICIT DOUBLE PRECISION! A~H, 0-2 ) 

Double precision XCK,M>,XC(NI 

0 — — — — 

DO 80 J « I , N 
i;,: XC( J ) = 0 , 0 

DO 10 IL =1 ,K1 
to XC!J)*XG(JF)+X!IL,J} 

II 1 ;; RK = K 1 

SO XC< J)*(XCtJ)-X!IEV1 , J> 1/CRK-1 .01 

; RETURN 

’ END 

C ************* ********m*******^****m*** *************** ****** 

SUBROUTINE CONSX ( N , H , K , ALPHA , BETA, GAMMA , DELTA , X , 

| 1 F , I EV£ » NO , G » H, XC* I PR I NT ) 

C COORDINATES SPECIAL PURPOSE SUBROUTINES 

1 : 6 ' 

C ARGUMENT LIST 

$€■"' 

c IT* ITERATION INDEX 

C IEV1 “INDEX OF POINT WITH HIN. FUNCTION VALUE 

C IEV£*IN0EX OF POINT WITH MAX. FUNCTION VALUE 

C I “POINT INDEX 

C KOBE* CONTROL KEY USED TO DETERMINE IF IMPLICIT CONSTRAITS 

C ARE PROVIDED 

C K 1 “DO LOOP LIMIT 

C ALL OTHHERS ARE PREVIOUSLY DEFINED IN MAIN LINE 

IMPLICIT DOUBLE PRECISION! A-H, G-Z I 
INTEGER GAMMA 

DIMENSION X!K,M>,F!KJ ,G!H>,H!HJ ,XC(NJ 
COMMON /B31/IT, ITHAX 
common /bb/i choice 
: ' IT*1 

f; , ■ 

I if ( i choice . eq.Dgo to 333 

DO 40 II *£,K 
DO 30 3*1 , N 
X!U,JJ*0.0 
38 CONTINUE 

48 CONTINUE 

C CALCULATE COMPtEX POINTS CHECK AGAINST CONSTRAINTS 

write! 6, *>' CALCULATING COMPLEX POINTS CHECK AGAINST CONSTRAINTS 
write{66 # *y*CALCULATTNS COMPLEX POINTS * 

&"l: ■ ■ do ss ''ti»»E'»Kl'V| ; ' 

rV-/-"' DO so J*!#N 

if ' • I»ll 'v.''.: ", ■ ' " ' , .,!■ 7 ■. ■. 

CALL CONST ! N , H , K , X, S ,H , 1 1 ;1 \ 

^ .;'r€0NT : iMUE^V7cl/ : ;^ 7 : r 



'$■ K 1 — 1 1 

“ < I . -£ ) " I S i N ss ' X ' 6 ' " • 1 ' K0DE ' XC ' DELTA ' K " 

51 IF( IPRINT) 52,65,52 

52 WRITECNO,018) 

WRITEC60, 018) 

ft 8 FORMAT ( // , 2X, 'CORDINATES OF 

10=1 

WRITE (NO, 01 9) (10, J,XCIO,J>, J=) ,N) 

WRITE(60 l 010MIO,J J ,XCI0 # J),Jaf,N) 

f!9 FORMAT ( / , 3 ( EX , 2HX C , 12, 1 H, , IE,4H) 

IFC IPRINT >56, 65,56 

56 WRITECNO, 019HII, J.XCII.J), j*l ,N> 

55 CONTINUE 

333 K 1 =>K 

DO 70 I *1 , K 
CALL FUNC (N,M,K,X,F, I) 

W CONTINUE 

KOUNT * 1 
IA*0 


INITIAL COMPLEX' ) 


* ,£ 13 . 6 )) 


FIND POINT OF LOWEST FUNCTION VALUE 
writeCfi,*) 'FIND POINT OF LOWEST FUNCTION VALUE' 
writ©(60,*> 'FIND POINT OF LOWEST FUNCTION VALUE' 


IF( IPRINT )7S, 80, ?£ 

72 WRITECNO, Q£1 ) 

WRITEC 60 , OST ) 

WRITEC 6 , 021 ) 

021 FORMAT C / , 2X , 22HVALUES OF THE FUNCTION ) 

WRITECNO, 022) CJ,PtJ),J*l ,K) 

WRITE C 60, OES) CJ,F(J) ,J*t , K > 

WRITEC 6 , 022JC J,FCJ), «I=t ,K> 

02S FORMAT C / , 3C EX , EHF C , 12 , 4H ) * ,1PEt3.6)J 

85 IEV1*1 

DO 1 00 I CM «e,K 
IFCFCIEV! )-F (I CM)) TO®, 100, SO 
35 IEV1=ICM 

150 CONTINUE . 

C FIND POINT WITH HI6HE8T FUNCTION VALUE 

WrlteCS,*) 'FIND POINT WITH HIGHEST FUNCTION VALUE' 
writeCSO, *) 'FIND POINT WITH HIGHEST FUNCTION VALUE* 

IEV2-1 

DO ISO I CM *2,K ";%&'& <f 

IFCFC IEVS)->FCieM) )f 10,110*12# 

IEV 2 *ICM :v r ■ 

120 .' CONTINUE ■ . V 

fc CHECK CONVERGENCE' CRITERIA 

wMtNCfi^ITCHECK COMVERSENCt CRITERIA* 
writ® (60,* ) 'CHECK CONVERGENCE CRITERIA' 



IF(F< IEV2 )-( F( IEV1 H-BETA) >140,130* 130 
130 KOUNT*1 

60 TO ISO 

140 KQUNT*KOUNT+1 

IF( K0UNT-6AHMA ) 150,240*240 

C REPLACE POINT WITH LOWEST FUNCTION VALUE 

writeCS,*)' REPLACE POINT WITH LOWEST FUNCTION VALUE' 
write? 60, * I 'REPLACE POINT WITH LOWEST FUNCTION VALUE' 

ISO CALL C£NTR(N*M,K, IEV1 , I,XC*X#K1 ) 

DO 160 JJ «=| , N 

160 X? IEV1 ,JJ)a t 1 . 0+ ALP HA )*?XC(JJ) 1— ALPHA* ( Xf IEV1 , JJ> ) 

161 I - I EV 1 

CALL CHECK ?N,H,K,X,G,H, I, KOBE, XC, DELTA* K1 > 

CALL FUNCCN,ri*K,X,F, I ) 

wr ite ? 6 , » ) ' F VALUE: ' , FCIEVti 

C REPLACE NEW POINT IF IT REPEATS AS LOWEST FUNCTION VALUE 

write?6,*) 'REPLACE NEW POINT REPEAT LOWEST FUNC VALUE' 
write? 60, *) 'REPLACE NEW POINT REPEAT LOWEST FUNC VALUE' 

170 IEV2*1 

DO ISO I CM *2,K 

IFtF? I£V2>-F(ICH) M90, 190, 180 
180 IEV2*ICH 

190 CONTINUE 

IF? IEV2-IEV1 >££0,200,220 
200 DO £100 JJ *1 ,N 

X? IEV1 , JJ)»?X( IEV1 , JJl+XC? JJ> >/2. 0 
2100 CONTINUE 

WRITE C 60 , *) 'CORDINATES OF ALL THE POINTS & F VALUES' 
do 9901 »r ® t ,K 

write? SO, *MXCi»r, Jcl , jc*t ,N), 'FC * *«r, ' I ' *FC»rl 
9901 continue 

162 X-1EV1 

CALL CHECKIN, H,K,X,C,H, I, KODE,XC, DELTA, K1 I 
CALL FUNCCN,«»k,X*F* I } 
write<6,*> 'F VALUE: ' , FC IEV1 > 

60 TO 170 
220 CONTINUE 

IFUPRINT J230, 228,230 
230 WRITE (NO, 023) IT 

WRITE? 60 , 023 ) IT 

023 FORMAT?//, 2X, 17HITERATI0N NUMBER ,I5> 

write? *,*> 'ItNo*' » IT 

WRITE? NO, 024 > 

WRITEC60, 024) 

024 FORMAT? //,2X, 'CORDINATES OF CORRECTED POINT') 

WRITE C NO , 01 9) f IEV1 , JC,X( lEVl , JC) * JC*1 *N) 

WRITE (60, 619 H IEVT * JC, X? IEV1 , JC ) , JC=1 ,N) 

WRI TE 1 6 0 i • 3> 'CORD I NATES OF ALL THE POINTS 4 F VALUES ' 

■ do 901 »r * 1, K ■ . r 



write(60, *> (X(mr, jc) , jc*t ,N) , 'Ft * ,ar, ' ) ' ,F(i*r) 

981 continue 

WRITE1NO, 021 ) 

WRITECNO,0£2)CI.F(I), I«1,K> 

WR I TE ( NO , OSS ) 

WR I TE ( 60 4 OS 1 ) 

WRITE (60, 02S ) 

925 FORMAT ( / # 2X , ' CORDINATES OF THE CENTROID' > 

WRITE (NO, 026) ( JC, XCC JC) , JC*1 , N) 

926 FORMAT ( / , 3 ( 2X , 2HX ( , I2,6H,C) « , 1 PEI 4.6, AX ) ) 

228 I T* I T+t 

IF ( I T- 1 THAX >89,80, 248 

240 RETURN 

END 

C ***************************************** 

SUBROUTINE CONST { N , M , K , X , C , H, I ) 

IMPLICIT DOUBLE PRECI SI ON ( A-H , 0-2 > 

Double precision 6(H) ,H(H),X(K,M) 

6(1 )» O.ldO 
H ( 1 >“ 560.0 
6(2)“ 1 .0 
H(2>» 5000.8 

6(3)“ 1.0 
H( 3 )*» 15000.0 
G<4)“ O.t 
H < 4 ) ® 10000.0 
6(5)“ 1.0 
H ( 5 )“ 45000.0 
6(6)* i .0 
H< 6 )“ 20000.0 
, G(7)“ 1 . 0 

H ( 7 )“ 50.0 
RETURN 
END 

c**m****m ******* ************************************************** 
SUBROUTINE FOKCC N, M, K * X, F, I ) 

IMPLICIT DOUBLE PRECISION* A-H, 0-2 ) 

DIMENSION X(K,M) ,F(K> 

COMMON /B4/ZETA1 ZETA2, ZETA3 
COMMON /B20/THETA1 * THETA2 , THETA3 
COMMON /8S/BW* TEMP, WS, DR 
COMMON /830/ERROR 

THETA1 * X( 1,11*1 . Od-4 
THETA2 ■ X( I , 2) 

THETA3 “ X(X,3) 

ZETA1 * X( I > 4>*t , Od-2 
ZETA2 “ XCI.S) 

■ ZETA3 « XC 1 , 6 ) ' .■ 1 

. Dll /. * X(I,?**t i0d“7 

. ' CALL SIMULATION - 



F C 2 > * ERROR 
write! * , * ) ' F VALUE* 
writet60, *) 'F VALUE 


ERROR 

ERROR 


RETURN 

END 

C** ******** ******** ******************* *m* *********9*** ttm****#*****i 
SUBROUTINE SIMULATION 
c using dQSebfe.f routine 

^s3sassasasBS*BesaaBSsasssBB3ssaBs=seB3sssasBBa*SBaaoaB»aa»ssass*B! 
C considering first reaction 
C fk£,rk£,fk3 are diffusion limited. 

C solvent has been considered. 

{JrasssasBBaaBBSSBBSSBissssBSBSBSBBissBaBBsaBBssrasraaaBSBasarBaaaaBa: 
IMPLICIT DOUBLE PRECISION (A-H,0-Z> 

DOUBLE PRECISION NINDEX 

* . . Parameters . . 

PARAMETER (N*66, IU*5500,NP«t 1 » 

PARAMETER (N0UT*6) 

* . . Local Arrays . . 

DOUBLE PRECISION Wt IW) , YIN) 

* . . External Functions . . 

EXTERNAL D0£EJW, S 

* .. External Subroutines 

EXTERNAL DOSEBF , D8EEJY, FCN, OUT, PEDER V , DOSCBF 

* .. Intrinsic Functions .. 

INTRINSIC DBLE 

c DIMENSION STATEMENTS 

c DIMENSION A0(3),E0(3>, ACC 3 ) » EC ( 3 ) , DH(3),BS(3> 

c COMMON STATEMENTS 

COMMON XEND, H, I 

COMMON /SE/R(Np) ,FKS£(Kp> ,FKS3(«p) ,RKS£INp>,RKS3CNp> 

COMMON /03/EKE, EK3 
COMMON /B4/ZETA t,ZETA£, ZETA3 
COMMON /BEO/THETAJ , THETAE, THETA3 
COMMON /B5/0II, TEMP , US, DR 
COMMON /B8/RR , P I 

COMMON /B9/THETA, AT, BT.TCINFIHITY, PI ,ALPHAL, ALPHAS 

COMMON /Bit /RAD I US 

COMMON /Bt 3/TS < NP ) , VFC NP ) , DDOCNP ) 

COMMON /B 1 4/FX t C NP ) , FX£ (HP), FX3 (NP), FX4 f NP l , FX5CNP ) , FX6CNP I 

COMMON /B1 5/AVEC1 , AVEP1 , AVELAMO, AVELAMt , AVELAM2, AVEW 

COMMON /Bi 6/AVEVF, AVETC 

COMMON /BI 8/VHS, VPS, GAMMA, VMO , ZETA, VFM 

COMMON /Bt 9 /NINDEX, WH, USGL 

COMMON /BEt/FKSHNP) ,RKSt CNPl 

COMMON /BEE/EKt 

COMMON /BE3/DD0E ( np ) , DD03 C np ) 

COMMON /B£4/VS0,VFS,VSS 
COMMON /BES/ZEtft 13, ZETAE3 

..■■v' COMNON /Bft/VlNIC0,6> ■Vv:. 



COMMON /B28/EXMUN ( 5 , 13) ,NS£T, I SET 
COMMON /B29/ERR0RINI , ICOUNT 
COMMON /B30/ERRQR, TOL 
COMMON /B31 /IT, UMAX 

ERROR * 0.000 
00 777 ISET ® 1 , N8ET 
ERROR INI = 0. OdO 


• INITIAL CONDITIONS .... 

DO JJ « 1 , Np 

* YINIC ISET, 1 ) 

Y(Np+JJ) a YINIC ISET , 2 ) 
Y<2*NP+JJ) = YINI ( ISET, 3) 
YC3*NP+JJ> - YINIC ISET, 4) 

Y ( 4*NP+ J J ) » YINIC ISET, S) 

Y(5*NP+JJ ) » YINIC ISET, 6 ) 

IFC JJ . EG . NP ) THEN 

Y C 5*NP+ JJ ) « 0 . 00001 dO 

ENDIF 

END DO 

writers,*) 'SIMULATION CALLED' 
writeCSO,*) 'SIMULATION CALLED' 

TEMP ■= 463 . 1 6D0 
TGINFINITY » 350 . ODO 
PI « O. 03DO 
ALPHAL ■ 1 . 0D-3 
ALPHAS « 4 .SD-4 
RADIUS « O. ©EDO 
MS * 0 . OOOOIdO 

C A_K PARAMETERS & CONSTANTS 

NINDEX * 8. OdO 
VMS ® 0.864E6300 
VPS * O.81232D0 
VSS » 0.969D0 
2ETA * 0 . 70924 6d0 
GAMMA a 1 . 000 
VM0 « 1 . 13247D0 
VFM » 0.31O33DO 
VS0 a t . 1429D0 
VFS * 0. 1795d0 

MM « t f 3. f 6D0 
MSOL a 18. ODO 
ZETA13 « 0 . 7O9296D0 
ZETA23 a 0.126500 

. X a 0,6B0 

EKt a a .35 IS 03 : .'V/..-:- 
c EKE * DEXPC C DSC EV-DHCS) / TEMP ) / R R ) 

C EK3 * BEXP ( C 08 <3 > -0H C 3 ) / TEMP ) /R R ) 

EKE a t 031.2474^0 

"• ■ ' . EiC3: ; :'a :-2, 422^0 ■' 



» 

DR * 0 . OOSdO 
DO JK » 1 , NP 
Rt JK) * DR* ( JK - f ) 

END DO 
■* 

* ■- Executable Statements .. 

c WRITE ( NOUT , * ) ' NYLON 6 SSP SIMULATION ' 

XEND *24 . ODO 
HPED * 0 
IR * 2 

I a 11 

H - ( XEND-X ) /DSLE (1+1 ) 

IFAIL « 0 
ICOUNT » 0 

CALL 002EBF ( X , XEND , N , Y , TOL, I R > FCN, HPEQ>D02£JY, OUT , U, III, IFAIL 
C CALL D02CBF ( X , XEND , N, Y , TOL , IR, FCN, OUT, W, IFAIL} 

IF t TOL . LT . 0 , ODO 1 WRITE (MOOT,*) ' Range too short for TOL * 

ERROR » ERROR + ERRORIN1 
777 CONTINUE 

55 RETURN 

writ e ( 6 , * ) ' F VALUE® * , ERROR 

99999 FORMAT <tX,A,D8.t> 

99998 FORMAT < IX, A, FT. 3) 

99997 FORMAT C 1 X , A , 6 C C 1 1 El 6 . 5 ) / ) ) 

99996 FORMAT < F4 . S/6C C 1 1 El 6 . 5 >/ ) ) 

ENO 

* — — 

SUBROUTINE FCNCT,Y,F) 

IMPLICIT DOUBLE PRECISION <A-H,0-Z) 

DOUBLE PRECISION N INDEX 
PARAMETER CN«66,Np»l 1 ) 

DIMENSION FCN), YCM) 

COMMON /B2/R ( Np ) , FKS2 C Np ) , FKS3C Np ) , RKS2 C Np > , RKS3C Np) 

COMMON /B5/DW , TEMP , US , OR 
COMMON /BE1/FKS1 (NP),RKSt (NP) 
common /const/const 
c . . Executable Statements . . 

DR a 0. OOSDO 

DO SS 33 » t , Np 

IJt « 33 

IJS \ •' a 1»Np + JJ 

133 , a 2*»p+JJ 

.. 134 ' . ' ■» 3*Np + 33 \ 

135 ' V' « ‘ '■ ■ 

, : ... : ; 

Cj^LL RATE CONST ( JJ , T : J3, 1 34 ,IJ6,Y> 



35 

25 


ZZZ * 
XL3 

n 

FC I Jt ) 


F t I J2 ) 


* 

* 

* 


Ft I J3) 


+FKS1 ( JJ )*Y( IJ1 )*Y< IJ6) - RKS1 ( j j)*YCIJ2J 
- Y 1 1 J5 ) * t £.0D0*Y(IJ5)*Y(IJ3) - Y( I J4)*Yf IJ4) 
/ < Y( I J4 ) * Y t I J3 ) ) 

“ -FK83t JJ)*Y( IJ1 ) * Y ( I J3 ) 

+RKS3tJJ)*tYt IJ3) - Yt IJE) ) - ZZZ 
* -2. 0DC*FKS2CJJ)*YC IJ2>*Y( IJ3i 

+2. 0DQ*RKS2< JJ>*Yt I J6)*(Yt I J3> - YUJfin 
~FKS3( JJ)*Y( IJ2)*YUJ1 ) 

+RK83 tJJ)*const*Yt IJ2> + ZZZ 
« -FKS2 tJJ)*Y(IJ3)*Y(IJ3) 


* +RKS2C JJJ*Y( IJ6>*tTt IJ4J - Yt I J3> ) 

* +ZZZ 

Ft IJ4) - FK83( JJ)*Yt IJ1 >*Yt IJ3) 

* -RK83 t J J ) * t Y t IJ3> - YCIJE)) 

* +ZZZ 

Ft IJ5) » 2. ODO*FKS£t JJ)*Yt IJ4)*Yt IJ4> 

* +1 .000/3. 0D0*RK82t JJ)*Yt IJ6)*tYt IJ4) -XL3I 

* +FKS3< JJ)*Yt IJ1 > * ( Yt I J3 ) + 2. 0D0*Yt IJ4) ) 

* +RKS3 ( JJ ) * t Y t IJ3i - 2. 0D0*Yf IJ4) + Yf IJ2I i 

* +ZZZ 

IFt JJ .EQ.NpJGO TO 35 
IFt JJ.EG. 1 ) THEN 
YP * Yt 57 ) 

Ft IJ6) * -Ft IJ3) 


* +Dw*ttYt 1J6+1 )-2.0D0*Yt IJ6M-YP>/tDR*BR) > 
ELSE 

YP » Y t i J 6-1 i 
F 1 1 J6 ) * -Ft IJ3) 

* +Dw*t f Y(I J6+I )-2. OB0*YUJ6I+YP>/t»R*DR J I 

* +D«*( Yt I J6+I ) - YP )/Rt JJ >/DR 
endi f 


GO TO 25 
YP * Y t 65 ) 
FCIJ6) * O.OB0 
CONTINUE 


» 



RETURN 

END 

C*aat**=*«»atss xaessss ssaastBss at seas as ass as ssass a»atiaaast5!at»ata*eaaassaasaiaa8aa*a!*!a«Bsx!B!aet>B«»aa»!B*ass 

SUBROUTINE OUTtX,Y) 

IMPLICIT DOUBLE PRECISION f A-H.O-Z) 

DOUBLE PRECISION N INDEX 
* . . Parameters 

PARAMETER (N«66,NP*t t ) 

PARAMETER tN0UT*6J 

*' ■ ■ . . Array Arguments .. ■. 

DOUBLE PRECISION YCN) 

'■ DOUBLE PR6.CI SION RRR^NF I . ■■ .. ■ 

tV^r , common ■ XEND* H , . ■ ■' • ■; . ■ - • ■. ■■ 

COMMON ■■ Mp I Up *VFMS3tNp *>RKS2(Np T>RK&3(M|» I 

^ : v ■ COMMON /BO-ZRADIUb;: ■>: ;.'V!r, ,• : ■■ V . , - . ■ . 

COMMON 



eOHHON /B14/FX1 C NP ) , FX2 ( NP ) , FX3 < NP ) , FX4 ( NP ) , FX5C NP ) , FX6 C HP > 
COMMON /B15/AVEC) , AVEP1 , AVELAMO , AVELAM1 , AVELAH2, AVEW 
COMMON /B16/AVEVF, AVETC 
COMMON /B17/FX7CNP) ,FX8(NP ) .FX9CNP) 

COMMON /B21/FKS1 (NP) ,RKS1 (NP> 

COMMON /B28/EXHUNC 5,13), N8ET , ISET 
COMMON /B29/ERRORINI, ICOUNT 
COMMON /B31/IT, ITHAX 


ICOUNT 

as 

ICOUNT 

+ 

DO JJ 

« 1 

.NP 



I Jt 

SB 

JJ 



IJS 

S3 

1 *Np 

+ 

JJ 

IJ3 

as 

2*Np 

+ 

JJ 

IJ4 

a 

3*Np 

* 

JJ 

IJ5 

a 

4*Np 

♦ 

JJ 

I J6 

a 

5*Np 


JJ 


RRR ( JJ ) » R ( JJ } *R ( J J ) 

CALL RATE CONST ( JJ, I J3, I J4 , I J6 , Y ) 

FXf < JJ) • RRR ( JJ ) *Y C I J1 ) 

FX2CJJ) » RRR { JJ ) * Y ( IJ2) 

FX3C JJ) « RRRC JJ)*Y( IJ3) 

FX4CJJ) » RRR ( JJ ) *Y( IJ4) 

FXSCJJ) - RRRC JJ)*YC IJS) 

FX6< JJ) * RRRC JJ ) *Y ( I J6 ) 

$ FX7CJJ) «■ AMUNCJJ) 

FXSCJJ) * RRRC JJ)*VFC JJ) 

| FXSCJJ) « RRRC JJ)*TCC JJ) 

% END DO 

CALL AVERASE C Y , AVEMUN ) 

WRITEC 6, * ) 'TIME® ' ,X, ' AVEMUN : '.AVEMUN 

IFCIT.6T. C ITRAX-I 8 ) ) WRITEC21 ,22)X, AVEMUN 

C WRITE C 22 . 23 )X . AVEC) , AVEP 1 » AVELAMO , AVELAMt . AVELAH2, AVEW 

22 FORMAT (4X,F10.3,4X,FI0.4) 

23 FORMAT ( 4X , F5 . 2 , EX, 6 ( Ft 3 . 8 ) ) 

c — — Calculating value of objective function — — 

FUN a -ABSC (AVEMUN - EXHUN ( I SET . I COUNT ) ) ) /AVEMUN 
ERROR INI * ERROR INI + FUN 

. c — ' — — — — * — ~ * — — • — — — , . 

X » XEND - DBLEM ) *H 

I « I - 1 

RETURN 

* ■ .. . ; ■ >' 

*9999 FORMAT CF10.4/6( C6Ft 6 . 8/SF! 6 . 8)/)) 

)\ . end ;• • ^ 

ijsv¥i 

^ ia?m,o-z) 



PARAMETER <N«66,Np*1 1 ) 

DIMENSION Y ( N ) , FKS20 ( NP ) , FKS30 ( NP ) , RKS20 ( NP ) 

COMMON /B2/R ( Np ) , FKS2 { Np ) , FKS3 ( Np ) * RKS2C Mp i , RKS3I Np J 

COMMON /B3/EK2 , EK3 

COMMON /B4/ZETA1 , 2ET A2 , ZETA3 

COMMON /B9/THETA , AT , BT , TG INF INI TY , P 1 , ALPHAL , ALPHA© 
COMMON /B 1 3/TG ( NP ) , VF ( NP > ,DDO(NP > 

COMMON /B20/THETA1 , THETA2 , THETA3 
COMMON /B21/FKS1 (NP) ,RKS1 (NP) 

COMMON /B22/EK1 

COMMON /B23/DB0£<np> ,DD03<np) 

COMMON /B25/ZETA1 3, ZETA23 

CALL DIFFUSIVITYMY, JJ, i j3,i|4, i j6) 

XYZ1 » THETAt /BD0( J J > 

XYZ2 * THETAE/DD02 ( J J > 

XYZ3 » THETA3/DDG3I JJ) 


FKS1CJJ) * ZETA1 ♦Y ( I J3) 
FK520 I J J ) « ZETAE*Y ( I J3 ) 
FKS30CJJ) « ZETA3*YC IJ3) 


RKS 1 ( JJ) ■ FK31 t JJ)/£K1 
RKS20 C J J ) * FKS20 < J J ) /EK2 

RKS3( JJ) » FKS30( JJJ/EK3 

C diffusion affected rate constants 

FKSE(JJ) “ t . ODO/U . 0DO/CFKS20( JJ) ) + XYZ1*Y( IJ3) > 
RKS2CJJ) * 1 . 9D0/I 1 . 0D0/CRKB2OC JJ) ) + XYZ2*YIIJ3)I 
FKS3CJJ) * 1 -0D0/C1 . 0D0/tFKS30C JJ) ) + XYZ3*YCIJ3)) 

RETURN 

END 


tssssaaB«sa*aa» 5 »s» 5 "« Bass3StSS!!SsaS!S!5!i!::5 ® sBs:s3!!SS:ss * s ® SSSSl 
SUBROUTINE BIFFUS1 VITY1 C Y , j j , i j 3, I j4 , i J6 ) 

IMPLICIT DOUBLE PRECISION (A-H,0*-Z> 

DOUBLE PRECISION NINDEX 
PARAMETER C M«66 » NP=* t! > 


DIMENSION Y(N),B0{NP),BB1 (NP>,DD2tNP> 

COMMON /BS/DW , TEMP , US , DR 

COMMON /B9/THETA, AT, BT,TGINFINITY, PI, ALPHAL, ALPHAG 

COMMON /BH/RADIUS 

COMMON /B13/TS t NP ),VFtNP>, DDO t NP ) , AHON t NP > 

COMMON /B 1 8/VMS , VPS, GAMMA , VM0 , ZET A , VFM 
COMMON /Bt 9/NI NDEX , UM , WSOL 
COMMON /B20/THETA1 
COMMON /B23/0D0£(np ) , DD03(np) 

COMMON /B24/VS0* VFS,VSS 
COMMON /Ba5/ZR**!3,ZETAe3 


AMUNtJJ) *YtJJ4)/YUJ3> ■. - ■ 

Tfei JJJ » I.eDO/n .ftOe/Te INFINITY + P1/AMUNC JJ)) 

VF* JJ) » O.OESd 0 + f ALP HAL -» ALPHAS)*! TEMP - TG( JJ) ) 

Bmuji bt*vfc jj) )> 







VP6 * VPS* ( 1 . OdO+VF { J J ) ) 

TV * CYC JJ)*UH*VMO + Y( I J3 ) *AHUN ( JJ > *WM*VPO + Y C IJ6 )*WSGL*VSOI 

PHIH = Y ( J J ) *WM*VHQ/ TV 

PHIS «* YC I J6>*USOL*VSO/TV 

PHIP * Y( IJ3)*AHUN! JJ)*UM*VP0/TV 

DM * 1 . OdQ/VHO 

BP * 1 . OdO/VPO 

0S » 1 . OdO/VSO 

VFP *> VFtJJ) 

XNUM * GAMMA* (DM*PHIM*VM8/ZETAt 3 + BS*PHI B*VSS/ZETA23 + 

* DP*PHIP*VPS ) 

XDEN « DM*PHIM*VMS*VFH + DS*PHI S*VSS*VFS +DP*PHIF*VPS*VFP 
CHI * XNUM/XDEN 

CHIO ■ CAHMA/VFP 

DDfJJ) * DEXP ( -CHI + CHIO) 

DD02C JJ) « DEXP ( ZETA23* ( -CHI+CHI 0 ) ) 

DD03C JJ) « BEXP ( ZETA1 3* ! "•CHI+CHI 0 ) ) 

DDKJJ) * AHUN! J J ) ** ( -N INDEX ) * DD( JJ) 


Qaiaso:: 


DBO( JJ)«DD1 ( JJ) 

RETURN 

END 

sa«ssasess«*sasies!ess»ses::ss,:*sss*B®«iBs::B:Bss*:xe:sisc»:Bssn»ssjs5ax»:ssB**B8SBasBsss«sessss»ssm*tt«*! 

SUBROUTINE AVERAGE ( Y , AVEMUN ) 

IMPLICIT DOUBLE PRECISION! A-H, O-Z ) 

DOUBLE PRECISION NINDEX 
PARAMETER (N«OB,NP=1t) 

DIMENSION YIN) 

COMMON /BT 1 /RADIUS 

COMMON /B13/TGINP) ,VF!NP> ,DB0!NP), AMUNtNP) 

COMMON /B 1 4/FX 1 ! NP ) , FX2 ! NP ) , FX3 t NP ) , FX4 f MP ) , FXS CUP ) , FX6! NP I 
COMMON /B 1 7/FX7 C NP ) » FXS (HP), FXS ( NP ) 

COMMON /B t S/AVEC1 , AVEP I , AVELAMO, AVELAMt , AVELAME , AVEM 
COMMON /B1 6/AVEVF, AVET6 


DENG » C 0 . 02tf 0**3 . OdO ) /3 . 0D0 
CALL INTEGRATION! FXt ,0.0,0.02, 1 1 , AINTt ) 
CALL INTEGRATION! FX2, 0.0,0.02,1 1 , AINT2) 
CALL INTEGRATISN!FX3, 0.0,0.02,1 t , AINT3) 
CALL INTEGRATION !FX4,0 .0,0.02,1 t , AINT4) 
CALL I NTEGR AT IONIFX5, O . 0, 0. 02, 1 1 , A I NTS ) 
CALL INTEGRATION! FX6 ,0.0,0.02,11 »AINT6) 
CALL INTEGRATI0N(FX7 ,0.0,0.02,11 , AINT7 ) 
CALL INTEGRATION! FXS, 0 .0,0.02,1 1 , AINTS ) 
CALL INTEGRATION! FX9, 0, 0,0. 02,1 1 , A I NTS J 
AVEC1 «AINTt/DE»0 
AVEP ! « AINT2/DEN0 

AVELAHO - AINT3/BEN0 
A VEL AM 1 * AI NT4/0ENG 



AVELAMS » AINT5/DENG 

J ^EW ^ » AINT6/DENQ 

WEVG I l § AINT8/DEN0 
VETS » AINT9/DENQ 
AVEMUN « AINT4/AINT3 

RETURN 

END 

SUBROUTINE INTEGRAT I QN ( FX , X I , XL, NPOINTS, XI NT) 

THIS SUBROUTINE INTEGRATES FUNCTION 'FX' OVER THE RANGE 
’ NPOINTS-1 ' NO. OF INTERVALS. NPOINTS SHOUO BE ODD! 
'XINT' IS THE VALUE OF INTEGRAL USING SIHPSONS' RULE. 

IMPLICIT DOUBLE PRECISION (A-H,0-2) 

PARAMETER (NP*t 1 ) 

DIMENSION FX(NP) 

N 1 » NPOINTS - 1 
NS » NPOINTS - S 
DX * 0.002 

SUM! » 0.0 
DO 3 J ,ss 3, NS, £ 

SUM1 * SUM! + FX(JJ) 

END DO 

SUMS * 0.0 

do jj « e,Ni,e 

SUMS « SUMS + FX(JJ) 

END DO 

SUM «* FXC ! ) + FX( NPOINTS) + S . 0D0*SUMf + 4. 0DO*SUH2 
XINT * DX*SU«/3.0 

RETURN 

, END; , 



